Bayesian GLMM Part2

Preparations

Load the necessary libraries

library(tidyverse)  #for data wrangling etc
library(rstanarm)   #for fitting models in STAN
library(cmdstanr)   #for cmdstan
library(brms)       #for fitting models in STAN
library(standist)   #for exploring distributions
library(HDInterval) #for HPD intervals
library(posterior)  #for posterior draws
library(coda)       #for diagnostics
library(bayesplot)  #for diagnostics
library(ggmcmc)     #for MCMC diagnostics
library(rstan)      #for interfacing with STAN
library(emmeans)    #for marginal means etc
library(broom)      #for tidying outputs
library(DHARMa)     #for residual diagnostics
library(tidybayes)  #for more tidying outputs
library(ggeffects)  #for partial plots
library(broom.mixed)#for tidying MCMC outputs
library(patchwork)  #for multiple plots
library(ggridges)   #for ridge plots 
library(bayestestR) #for ROPE
library(see)        #for some plots
source('helperFunctions.R')

Scenario

To investigate differential metabolic plasticity in barramundi (Lates calcarifer), Norin, Malte, and Clark (2015) exposed juvenile barramundi to various environmental changes (increased temperature, decreased salinity and increased hypoxia) as well as control conditions. Metabolic plasticity was calculated as the percentage difference in standard metabolic rate between the various treatment conditions and the standard metabolic rate under control conditions. They were interested in whether there was a relationship between metabolic plasticity and typical (control) metabolism and how the different treatment conditions impact on this relationship.

A total of 60 barramundi juveniles were subject to each of the three conditions (high temperature, low salinity and hypoxia) in addition to control conditions. Fish mass was also recorded as a covariate as this is known to influence metabolic parameters.

Barramundi

Format of norin.csv data files

FISHID	MASS	TRIAL	SMR_contr	CHANGE
1	35.69	LowSalinity	5.85	-31.92
2	33.84	LowSalinity	6.53	2.52
3	37.78	LowSalinity	5.66	-6.28
..	..	..	..	..
1	36.80	HighTemperature	5.85	18.32
2	34.98	HighTemperature	6.53	19.06
3	38.38	HighTemperature	5.66	19.03
..	..	..	..	..
1	45.06	Hypoxia	5.85	-18.61
2	43.51	Hypoxia	6.53	-5.37
3	45.11	Hypoxia	5.66	-13.95

FISHID	Categorical listing of the individual fish that are repeatedly sampled
MASS	Mass (g) of barramundi. Covariate in analysis
TRIAL	Categorical listing of the trial (LowSalinity: 10ppt salinity; HighTemperature: 35 degrees; Hypoxia: 45% air-sat. oxygen.
SMR_contr	Standard metabolic rate (mg/h/39.4 g of fish) under control trial conditions (35 ppt salinity, 29 degrees, normoxia)
CHANGE	Percentage difference in Standard metabolic rate (mg/h/39.4 g of fish) between Trial conditions and control adjusted for 'regression to the mean'.

Read in the data

norin <- read_csv('../public/data/norin.csv', trim_ws=TRUE)

## Rows: 180 Columns: 5
## ── Column specification ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): TRIAL
## dbl (4): FISHID, MASS, SMR_contr, CHANGE
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

glimpse(norin)

## Rows: 180
## Columns: 5
## $ FISHID    <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 1…
## $ MASS      <dbl> 35.69, 33.84, 37.78, 26.58, 37.62, 37.68, 30.62, 50.37, 24.9…
## $ TRIAL     <chr> "LowSalinity", "LowSalinity", "LowSalinity", "LowSalinity", …
## $ SMR_contr <dbl> 5.847466, 6.530707, 5.659556, 6.278200, 4.407336, 4.818589, …
## $ CHANGE    <dbl> -31.919389, 2.520929, -6.284968, -4.346675, -3.071329, -15.0…

We should start by declaring any categorical variables (including an varying effects) as predictors.

Exploratory data analysis

Model formula: \[ y_i \sim{} \mathcal{N}(\mu_i, \sigma^2)\\ \mu_i =\boldsymbol{\beta} \bf{X_i} + \boldsymbol{\gamma} \bf{Z_i}\\ \beta_0 \sim{} \mathcal{N}(16, 35)\\ \beta_{1-6} \sim{} \mathcal{N}(0, 70)\\ \]

where \(\boldsymbol{\beta}\) and \(\boldsymbol{\gamma}\) are vectors of the fixed and random effects parameters respectively and \(\bf{X}\) is the model matrix representing the overall intercept and effects of temperature and (centered) mean fish size on SDA peak. \(\bf{Z}\) represents a cell means model matrix for the random intercepts associated with individual fish.

Lets begin by declaring the categorical predictors and random effects as factors.

norin <- norin %>% mutate(FISHID=factor(FISHID),
                         TRIAL=factor(TRIAL))

Each of the fish were exposed to each of the three trials (treatments). So lets explore the distributions of responses within each of these trials.

ggplot(norin, aes(y=CHANGE, x=TRIAL)) + geom_boxplot()

Conclusions:

• these seem reasonable enough

The researchers considered that physiological plasticity might be effected by the fishes basal metabolic rate. For example, a fish with relatively high metabolism might have less scope to increase this metabolism than a fish with a lower metabolism. Therefore, the SMR under control conditions (just prior to the specific trial) could be added as a continuous covariate.

Doing so would introduce two additional assumptions:

1. linearity: that the relationship between the response and control SMR is linear (if we intend to model it as a linear trend)
2. homogeneity of slopes: that the rate of change in response associated with a one unit change in control SMR is similar for each trial (unless we include an interaction)

ggplot(norin, aes(y=CHANGE, x=SMR_contr, shape=TRIAL, color=TRIAL)) +
    geom_smooth(method='lm') + geom_point()

## `geom_smooth()` using formula 'y ~ x'

ggplot(norin, aes(y=CHANGE, x=SMR_contr, shape=TRIAL, color=TRIAL)) +
  geom_smooth() + geom_point()

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

It might also be worth exploring how consistent the trial effect is across the different fish. This can give us an idea of whether the addition of a random intercept/slope model would be useful.

ggplot(norin, aes(y=CHANGE, x=as.numeric(FISHID), color=TRIAL)) +
    geom_point() + geom_line()

Conclusions:

• generally, the difference between control and trial SMR is greatest for the High temperature trial.
• there is less consistency about the relative effects of Hypoxia and Low Salinity.
• hence a random intercept/slope model might be useful.

Finally, as an acknowledgement that the metabolic response might be influenced

by the mass of the fish, the researchers contemplated including fish Mass as a continuous covariate. There are numerous alternative ways to incorporate covariates that are known to impact on a response.

• add them as an offset. This will include the covariate in the model, yet will not estimate a partial slope for this predictor. Rather the slope is assumed to be 1. Since the parameter does not need to be estimated, adding a predictor as an offset does not consume any degrees of freedom. However, it does assume that there is a slope of one. This approach can be useful as an alternative to using the response divided by the covariate as the response. That is, it is equivalent to using a response that has been standardised by the covariate, yet still allowing the modelling distribution that would be appropriate for the original response.
• add them as a regular covariate. This will include the covariate in the model and estimate a partial slope parameter just like any other term. Although it will incur a degrees of freedom penalty, it does not assume a slope of 1 and could be modelled as non-linear if appropriate.

Lets explore the relationship between the response and Mass separately for each Trial.

ggplot(norin, aes(y=CHANGE, x=MASS, color=TRIAL)) +
  geom_point() +
  geom_smooth(method='lm')

## `geom_smooth()` using formula 'y ~ x'

Conclusions:

• there does not seem to be much of a relationship between the response and Mass. Nevertheless, in the presence of control SMN, it might still be a useful covariate at explaining some of the unexplained variability (and thus increasing the power of the model).
• there is no evidence that the response and Mass are related in a non-linear manner.

Fit the model

rstanarm

Using default priors

In rstanarm, the default priors are designed to be weakly informative. They are chosen to provide moderate regularlization (to help prevent overfitting) and help stabalise the computations.

norin.rstanarm <- stan_glmer(CHANGE ~  (1|FISHID)+TRIAL*scale(SMR_contr, scale = FALSE)+scale(MASS, scale = FALSE),
                             data = norin,
                             family = gaussian(), 
                             iter = 5000,
                             warmup = 2000,
                             chains = 3,
                             thin = 5,
                             refresh = 0) 

norin.rstanarm %>% prior_summary()

## Priors for model '.' 
## ------
## Intercept (after predictors centered)
##   Specified prior:
##     ~ normal(location = 20, scale = 2.5)
##   Adjusted prior:
##     ~ normal(location = 20, scale = 85)
## 
## Coefficients
##   Specified prior:
##     ~ normal(location = [0,0,0,...], scale = [2.5,2.5,2.5,...])
##   Adjusted prior:
##     ~ normal(location = [0,0,0,...], scale = [178.99,178.99,135.15,...])
## 
## Auxiliary (sigma)
##   Specified prior:
##     ~ exponential(rate = 1)
##   Adjusted prior:
##     ~ exponential(rate = 0.03)
## 
## Covariance
##  ~ decov(reg. = 1, conc. = 1, shape = 1, scale = 1)
## ------
## See help('prior_summary.stanreg') for more details

This tells us:

• for the intercept (when the family is Gaussian), a normal prior with a mean of 20 and a standard deviation of 85 is used. A value of 2.5 is used for scaling all parameter standard deviations. The value of 20 is based on the mean of the response (mean(norin$CHANGE)) and the scaled standard deviation of 85 is based on multiplying the scaling factor by the standard deviation of the response.

2.5*sd(norin$CHANGE)

## [1] 84.6135

• for the coefficients (in this case, just the difference between strong and weak innoculation), the default prior is a normal prior centered around 0 with a standard deviations of 178.99, 178.99, 135.15, 10.60, 34.2 and 34.2 repectively. This is then adjusted for the scale of the data by multiplying the 2.5 by the ratio of the standard deviation of the response by the stanard devation of the numerical dummy variables for the predictor (then rounded).

2.5*sd(norin$CHANGE)/apply(model.matrix(~TRIAL*SMR_contr+MASS, norin), 2, sd)

##                (Intercept)               TRIALHypoxia 
##                        Inf                  178.99305 
##           TRIALLowSalinity                  SMR_contr 
##                  178.99305                  135.14517 
##                       MASS     TRIALHypoxia:SMR_contr 
##                   10.59372                   34.19156 
## TRIALLowSalinity:SMR_contr 
##                   34.19156

• the auxillary prior represents whatever additional prior is required for the nominated model. In the case of a Gaussian model, this will be \(\sigma\), for negative binomial, it will be the reciprocal of dispersion, for gamma, it will be shape, etc . By default in rstanarm, this is a exponential with a rate of 1 which is then adjusted by devision with the standard deviation of the response.

1/sd(norin$CHANGE)

## [1] 0.02954611

Assessing priors

Lets now run with priors only so that we can explore the range of values they allow in the posteriors.

norin.rstanarm1 <- update(norin.rstanarm,  prior_PD=TRUE)

norin.rstanarm1 %>%
  ggpredict(~SMR_contr*TRIAL) %>%
  plot(add.data=TRUE)

## Note: uncertainty of error terms are not taken into account. You may want to use `rstantools::posterior_predict()`.

#OR
norin.rstanarm1 %>%
  ggemmeans(~SMR_contr*TRIAL) %>%
  plot(add.data=TRUE)

Conclusions:

• we see that the range of predictions is faily wide and the predicted means could range from a small negative number to a relatively large positive number.

Defining priors

The following link provides some guidance about defining priors. [https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations]

When defining our own priors, we typically do not want them to be scaled.

If we wanted to define our own priors that were less vague, yet still not likely to bias the outcomes, we could try the following priors (mainly plucket out of thin air):

• \(\beta_0\): normal centered at 17 with a standard deviation of 35
  • mean of 17: since median(norin$CHANGE)
  • sd of 35: since mad(norin$CHANGE)
• \(\beta_1\): normal centred at 0 with a standard deviation of 70
  • sd of 70: since sd(norin$CHANGE) / apply(model.matrix(~scale(SMR_contr)*TRIAL + scale(MASS), norin), 2, sd)
• \(\sigma\): exponential with rate of 0.03
  • since 1 / sd(norin$CHANGE)
• \(\Sigma\): decov with:
  • regularization: the exponent for a LKJ prior on the correlation matrix. A value of 1 (default) implies a joint uniform prior
  • concentration: the concentration parameter for a symmetric Dirichlet distribution. A value of 1 (default) implies a joint uniform distribution
  • shape and scale: the shape and scale parameters for a gamma prior on the scale and scale parameters of the decov prior. A value of 1 for both (default) simplifies the gamma prior to a unit-exponential distribution.

I will also overlay the raw data for comparison.

norin.rstanarm2 <- stan_glmer(CHANGE ~  (1|FISHID)+TRIAL*scale(SMR_contr, scale = FALSE)+offset(MASS),
                                data = norin,
                                family = gaussian(), 
                                prior_intercept = normal(17, 35, autoscale = FALSE),
                                prior = normal(0, 70, autoscale = FALSE),
                                prior_aux=rstanarm::exponential(0.03, autoscale = FALSE),
                                prior_covariance = decov(1, 1, 1, 1), 
                                prior_PD = TRUE, 
                                iter = 5000,
                                warmup = 1000,
                                chains = 3,
                                thin = 5,
                                refresh = 0
                                )  

norin.rstanarm2 %>%
    ggpredict(~SMR_contr * TRIAL) %>%
    plot(add.data = TRUE)

## Warning in sweep(eta, 2L, offset, `+`): STATS is longer than the extent of
## 'dim(x)[MARGIN]'

## Note: uncertainty of error terms are not taken into account. You may want to use `rstantools::posterior_predict()`.

Now lets refit, conditioning on the data.

norin.rstanarm3 <- update(norin.rstanarm2,  prior_PD=FALSE)

norin.rstanarm4 <- stan_glmer(CHANGE ~  (TRIAL|FISHID)+TRIAL*SMR_contr + offset(MASS),
                                data = norin,
                                family = gaussian(), 
                                prior_intercept = normal(17, 35, autoscale = FALSE),
                                prior = normal(0, 70, autoscale = FALSE),
                                prior_aux=rstanarm::exponential(0.03, autoscale = FALSE),
                                prior_covariance = decov(1, 1, 1, 1), 
                                iter = 5000,
                                warmup = 1000,
                                chains = 3,
                                thin = 5,
                                refresh = 0
                                )  

## Warning: There were 1 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.

## Warning: There were 3 chains where the estimated Bayesian Fraction of Missing Information was low. See
## https://mc-stan.org/misc/warnings.html#bfmi-low

## Warning: Examine the pairs() plot to diagnose sampling problems

preds <- norin.rstanarm4 %>% posterior_predict(ndraws = 250,  summary = FALSE)
norin.resids <- createDHARMa(simulatedResponse = t(preds),
                            observedResponse = norin$CHANGE,
                            fittedPredictedResponse = apply(preds, 2, median),
                            integerResponse = FALSE)
plot(norin.resids, quantreg = FALSE)

## Clearly an issue here!

Plotting prior and posterior

get_variables(norin.rstanarm3)

##  [1] "(Intercept)"                                     
##  [2] "TRIALHypoxia"                                    
##  [3] "TRIALLowSalinity"                                
##  [4] "scale(SMR_contr, scale = FALSE)"                 
##  [5] "TRIALHypoxia:scale(SMR_contr, scale = FALSE)"    
##  [6] "TRIALLowSalinity:scale(SMR_contr, scale = FALSE)"
##  [7] "b[(Intercept) FISHID:1]"                         
##  [8] "b[(Intercept) FISHID:2]"                         
##  [9] "b[(Intercept) FISHID:3]"                         
## [10] "b[(Intercept) FISHID:4]"                         
## [11] "b[(Intercept) FISHID:5]"                         
## [12] "b[(Intercept) FISHID:6]"                         
## [13] "b[(Intercept) FISHID:7]"                         
## [14] "b[(Intercept) FISHID:8]"                         
## [15] "b[(Intercept) FISHID:9]"                         
## [16] "b[(Intercept) FISHID:10]"                        
## [17] "b[(Intercept) FISHID:11]"                        
## [18] "b[(Intercept) FISHID:12]"                        
## [19] "b[(Intercept) FISHID:13]"                        
## [20] "b[(Intercept) FISHID:14]"                        
## [21] "b[(Intercept) FISHID:15]"                        
## [22] "b[(Intercept) FISHID:16]"                        
## [23] "b[(Intercept) FISHID:17]"                        
## [24] "b[(Intercept) FISHID:18]"                        
## [25] "b[(Intercept) FISHID:19]"                        
## [26] "b[(Intercept) FISHID:20]"                        
## [27] "b[(Intercept) FISHID:21]"                        
## [28] "b[(Intercept) FISHID:22]"                        
## [29] "b[(Intercept) FISHID:23]"                        
## [30] "b[(Intercept) FISHID:24]"                        
## [31] "b[(Intercept) FISHID:25]"                        
## [32] "b[(Intercept) FISHID:26]"                        
## [33] "b[(Intercept) FISHID:27]"                        
## [34] "b[(Intercept) FISHID:28]"                        
## [35] "b[(Intercept) FISHID:29]"                        
## [36] "b[(Intercept) FISHID:30]"                        
## [37] "b[(Intercept) FISHID:31]"                        
## [38] "b[(Intercept) FISHID:32]"                        
## [39] "b[(Intercept) FISHID:33]"                        
## [40] "b[(Intercept) FISHID:34]"                        
## [41] "b[(Intercept) FISHID:35]"                        
## [42] "b[(Intercept) FISHID:36]"                        
## [43] "b[(Intercept) FISHID:37]"                        
## [44] "b[(Intercept) FISHID:38]"                        
## [45] "b[(Intercept) FISHID:39]"                        
## [46] "b[(Intercept) FISHID:40]"                        
## [47] "b[(Intercept) FISHID:41]"                        
## [48] "b[(Intercept) FISHID:42]"                        
## [49] "b[(Intercept) FISHID:43]"                        
## [50] "b[(Intercept) FISHID:44]"                        
## [51] "b[(Intercept) FISHID:45]"                        
## [52] "b[(Intercept) FISHID:46]"                        
## [53] "b[(Intercept) FISHID:47]"                        
## [54] "b[(Intercept) FISHID:48]"                        
## [55] "b[(Intercept) FISHID:49]"                        
## [56] "b[(Intercept) FISHID:50]"                        
## [57] "b[(Intercept) FISHID:51]"                        
## [58] "b[(Intercept) FISHID:52]"                        
## [59] "b[(Intercept) FISHID:53]"                        
## [60] "b[(Intercept) FISHID:54]"                        
## [61] "b[(Intercept) FISHID:55]"                        
## [62] "b[(Intercept) FISHID:56]"                        
## [63] "b[(Intercept) FISHID:57]"                        
## [64] "b[(Intercept) FISHID:58]"                        
## [65] "b[(Intercept) FISHID:59]"                        
## [66] "b[(Intercept) FISHID:60]"                        
## [67] "sigma"                                           
## [68] "Sigma[FISHID:(Intercept),(Intercept)]"           
## [69] "accept_stat__"                                   
## [70] "stepsize__"                                      
## [71] "treedepth__"                                     
## [72] "n_leapfrog__"                                    
## [73] "divergent__"                                     
## [74] "energy__"

posterior_vs_prior(norin.rstanarm3, color_by='vs', group_by=TRUE,
                   facet_args=list(scales='free_y'),
                   regex_pars = "^.Intercept|TRIAL|SMR_contr|MASS|sigma|Sigma")

## 
## Drawing from prior...

Conclusions:

• in each case, the prior is substantially wider than the posterior, suggesting that the posterior is not biased towards the prior.

norin.rstanarm3 %>%
    ggpredict(~SMR_contr * TRIAL) %>%
    plot(add.data = TRUE)

## Warning in sweep(eta, 2L, offset, `+`): STATS is longer than the extent of
## 'dim(x)[MARGIN]'

## Note: uncertainty of error terms are not taken into account. You may want to use `rstantools::posterior_predict()`.

##OR
norin.rstanarm3 %>%
    ggemmeans(~SMR_contr * TRIAL) %>%
    plot(add.data = TRUE)

brms

Using default priors

In brms, the default priors are designed to be weakly informative. They are chosen to provide moderate regularlization (to help prevent overfitting) and help stabalise the computations.

Unlike rstanarm, brms models must be compiled before they start sampling. For most models, the compilation of the stan code takes around 45 seconds.

norin.form <- bf(CHANGE ~  (1|FISHID)+TRIAL*SMR_contr+offset(MASS),
                   family = gaussian() 
                   )
options(width=100)
norin.form %>% get_prior(data=norin)

##                                  prior     class                       coef  group resp dpar nlpar
##                                 (flat)         b                                                  
##                                 (flat)         b                  SMR_contr                       
##                                 (flat)         b               TRIALHypoxia                       
##                                 (flat)         b     TRIALHypoxia:SMR_contr                       
##                                 (flat)         b           TRIALLowSalinity                       
##                                 (flat)         b TRIALLowSalinity:SMR_contr                       
##  student_t(3, -23.7797222222222, 34.5) Intercept                                                  
##                  student_t(3, 0, 34.5)        sd                                                  
##                  student_t(3, 0, 34.5)        sd                            FISHID                
##                  student_t(3, 0, 34.5)        sd                  Intercept FISHID                
##                  student_t(3, 0, 34.5)     sigma                                                  
##  lb ub       source
##             default
##        (vectorized)
##        (vectorized)
##        (vectorized)
##        (vectorized)
##        (vectorized)
##             default
##   0         default
##   0    (vectorized)
##   0    (vectorized)
##   0         default

norin.brm <- brm(norin.form,
                 data=norin,
                 iter = 5000,
                 warmup = 1000,
                 chains = 3, cores = 3,
                 thin = 5,
                 refresh = 0,
                 backend = 'cmdstanr')

## Start sampling

## Running MCMC with 3 parallel chains...
## 
## Chain 2 finished in 3.3 seconds.
## Chain 3 finished in 3.2 seconds.
## Chain 1 finished in 3.5 seconds.
## 
## All 3 chains finished successfully.
## Mean chain execution time: 3.3 seconds.
## Total execution time: 3.7 seconds.

This tells us:

• for the intercept (when the family is Gaussian), a t-distribution prior with a mean of 16.4 and a standard deviation of 34.5 is used:
  • mean of 16.4, since median(norin$CHANGE)
  • sd of 34.5 since mad(norin$CHANGE)
• brms uses flat (inproper) priors for all population effects
• the prior on sigma is a half t-distribution with a mean of 0 and standard deviation of 34.5:
  • sd of 34.5 since mad(norin$CHANGE)
• the hyperprior on the standard deviation is a half t-distribution with a mean of 0 and a standard deviation of 34.5.

Defining priors

The following link provides some guidance about defining priors. [https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations]

When defining our own priors, we typically do not want them to be scaled.

If we wanted to define our own priors that were less vague, yet still not likely to bias the outcomes, we could try the following priors (mainly plucket out of thin air):

• \(\beta_0\): normal centred at 16 with a standard deviation of 35
  • mean of 16: since median(norin$CHANGE)
  • sd of 35: since mad(norin$CHANGE)
• \(\beta_{1-2}\) (TRIAL): normal centred at 0 with a standard deviation of 70
  • sd of 70: since sd(norin$CHANGE) / apply(model.matrix(~TRIAL*SMR_contr+MASS, norin)[, -1], 2, sd)
• \(\beta_{3}\) (SMR_contr): normal centred at 0 with a standard deviation of 54
  • sd of 54: since sd(norin$CHANGE) / apply(model.matrix(~TRIAL*SMR_contr+MASS, norin)[, -1], 2, sd)
• \(\beta_{4}\) (MASS): normal centred at 0 with a standard deviation of 4
  • sd of 4: since sd(norin$CHANGE) / apply(model.matrix(~TRIAL*SMR_contr+MASS, norin)[, -1], 2, sd)
• \(\beta_{5-6}\) (interaction): normal centred at 0 with a standard deviation of 14
  • sd of 14: since sd(norin$CHANGE) / apply(model.matrix(~TRIAL*SMR_contr+MASS, norin)[, -1], 2, sd)
• \(\sigma\): half-cauchy with parameters 0 and 6 or else a student_t with parameters 0 and 34
  • since sd(norin$CHANGE)
• \(\sigma_j\): half-cauchy with parameters 0 and 5.8
  • since sqrt(sd(norin$CHANGE))
  • we want this prior to have most mass close to zero for the purpose of regularisation
• \(\Sigma\): decov with:
  • regularization: the exponent for a LKJ prior on the correlation matrix. A value of 1 (default) implies a joint uniform prior
  • concentration: the concentration parameter for a symmetric Dirichlet distribution. A value of 1 (default) implies a joint uniform distribution
  • shape and scale: the shape and scale parameters for a gamma prior on the scale and scale parameters of the decov prior. A value of 1 for both (default) simplifies the gamma prior to a unit-exponential distribution.

Note, for hierarchical models, the model will tend to want to have a large \(sigma\) in order to fit the data better. It is a good idea to regularise this tendency by applying a prior that has most mass around zero. Suitable candidates include:

• half-t: as the degrees of freedom approach infinity, this will approach a half-normal
• half-cauchy: this is essentially a half-t with 1 degree of freedom
• exponential

I will also overlay the raw data for comparison.

norin %>% 
    group_by(TRIAL) %>%
    summarise(median(CHANGE),
              mad(CHANGE))

norin %>% 
    group_by(TRIAL, FISHID) %>%
    summarise(median = median(CHANGE),
              MAD = mad(CHANGE)) %>%
    ungroup(FISHID) %>%
    summarise(sd(median))

## `summarise()` has grouped output by 'TRIAL'. You can override using the
## `.groups` argument.

sd(norin$CHANGE)/apply(model.matrix(~TRIAL*SMR_contr+MASS, norin)[, -1], 2, sd)

##               TRIALHypoxia           TRIALLowSalinity 
##                  71.597218                  71.597218 
##                  SMR_contr                       MASS 
##                  54.058070                   4.237486 
##     TRIALHypoxia:SMR_contr TRIALLowSalinity:SMR_contr 
##                  13.676626                  13.676626

standist::visualize("normal(16,35)", xlim=c(-10,100))

standist::visualize("normal(0, 70)", "normal(0, 54)", xlim=c(-200,200))

standist::visualize("gamma(2, 1)", "gamma(35, 1)",
                    "student_t(3,0, 34)",
                    "cauchy(0, 5.8)",
                    xlim=c(-10,100))

priors <- prior(normal(16, 35), class = 'Intercept') +
    prior(normal(0, 70), class = 'b', coef = 'TRIALHypoxia') +
    prior(normal(0, 70), class = 'b', coef = 'TRIALLowSalinity') +
    prior(normal(0, 54), class = 'b', coef = 'SMR_contr') +
    prior(normal(0, 4), class = 'b', coef = 'MASS') +
    prior(normal(0, 13), class = 'b') +
    prior(student_t(3,0,34), class = 'sigma') +
    ## prior(gamma(35, 1), class = 'sigma') +
    prior(cauchy(0, 5.8), class = 'sd') 
norin.form <- bf(CHANGE ~ (1|FISHID) + TRIAL*SMR_contr + MASS,
                     family = gaussian()
                   )
norin.brm2 <- brm(norin.form, 
                  data = norin,
                  prior = priors,
                  sample_prior = 'only',
                  iter = 5000,
                  warmup = 1000,
                  chains = 3, cores = 3,
                  thin = 5,
                  refresh = 0,
                  backend = "cmdstanr"
                  )

## Start sampling

## Running MCMC with 3 parallel chains...
## 
## Chain 1 finished in 0.2 seconds.
## Chain 2 finished in 0.2 seconds.
## Chain 3 finished in 0.2 seconds.
## 
## All 3 chains finished successfully.
## Mean chain execution time: 0.2 seconds.
## Total execution time: 0.4 seconds.

Note in the above model, the output may have included a warning message alerting us the presence of divergent transitions. Divergent transitions are an indication that the sampler has encountered poor sampling conditions - the more divergent transitions, the more severe the issue.

Typically, divergent transitions are the result of either:

• a miss-specified model
• priors that permit the sampler to drift into unsupported areas
• complex posterior “features” (with high degrees of curvature) for which the sampler was inadequately tuned during the warmup phase

Accordingly, these divergent transitions can be addressed by either:

• reviewing the model structure
• adopting tighter priors
• increase the adaptive delta from the default of 0.8 to closer to 1. The adaptive delta defines the average acceptance probability that the sampler should aspire to during the warmup phase. Increasing the adaptive delta results in a smaller step size (and thus fewer divergences and more robust samples) however it will also result in slower sampling speeds.

norin.brm2 %>%
    ggpredict(~SMR_contr*TRIAL) %>%
    plot(add.data = TRUE)

## Note: uncertainty of error terms are not taken into account. You may want to use `rstantools::posterior_predict()`.

norin.brm3 <- update(norin.brm2,  
                     sample_prior = 'yes',
                     control = list(adapt_delta = 0.99),
                     refresh = 0,
                     cores = 3)  

## The desired updates require recompiling the model

## Start sampling

## Running MCMC with 3 parallel chains...
## 
## Chain 2 finished in 5.6 seconds.
## Chain 3 finished in 6.6 seconds.
## Chain 1 finished in 9.5 seconds.
## 
## All 3 chains finished successfully.
## Mean chain execution time: 7.2 seconds.
## Total execution time: 9.6 seconds.

save(norin.brm3, file = '../ws/testing/norin.brm3')

norin.brm3 %>%
    ggpredict(~SMR_contr*TRIAL) %>%
    plot(add.data = TRUE)

## Note: uncertainty of error terms are not taken into account. You may want to use `rstantools::posterior_predict()`.

Plotting prior and posterior

norin.brm3 %>% get_variables()

##  [1] "b_Intercept"                        "b_TRIALHypoxia"                    
##  [3] "b_TRIALLowSalinity"                 "b_SMR_contr"                       
##  [5] "b_MASS"                             "b_TRIALHypoxia:SMR_contr"          
##  [7] "b_TRIALLowSalinity:SMR_contr"       "sd_FISHID__Intercept"              
##  [9] "sigma"                              "r_FISHID[1,Intercept]"             
## [11] "r_FISHID[2,Intercept]"              "r_FISHID[3,Intercept]"             
## [13] "r_FISHID[4,Intercept]"              "r_FISHID[5,Intercept]"             
## [15] "r_FISHID[6,Intercept]"              "r_FISHID[7,Intercept]"             
## [17] "r_FISHID[8,Intercept]"              "r_FISHID[9,Intercept]"             
## [19] "r_FISHID[10,Intercept]"             "r_FISHID[11,Intercept]"            
## [21] "r_FISHID[12,Intercept]"             "r_FISHID[13,Intercept]"            
## [23] "r_FISHID[14,Intercept]"             "r_FISHID[15,Intercept]"            
## [25] "r_FISHID[16,Intercept]"             "r_FISHID[17,Intercept]"            
## [27] "r_FISHID[18,Intercept]"             "r_FISHID[19,Intercept]"            
## [29] "r_FISHID[20,Intercept]"             "r_FISHID[21,Intercept]"            
## [31] "r_FISHID[22,Intercept]"             "r_FISHID[23,Intercept]"            
## [33] "r_FISHID[24,Intercept]"             "r_FISHID[25,Intercept]"            
## [35] "r_FISHID[26,Intercept]"             "r_FISHID[27,Intercept]"            
## [37] "r_FISHID[28,Intercept]"             "r_FISHID[29,Intercept]"            
## [39] "r_FISHID[30,Intercept]"             "r_FISHID[31,Intercept]"            
## [41] "r_FISHID[32,Intercept]"             "r_FISHID[33,Intercept]"            
## [43] "r_FISHID[34,Intercept]"             "r_FISHID[35,Intercept]"            
## [45] "r_FISHID[36,Intercept]"             "r_FISHID[37,Intercept]"            
## [47] "r_FISHID[38,Intercept]"             "r_FISHID[39,Intercept]"            
## [49] "r_FISHID[40,Intercept]"             "r_FISHID[41,Intercept]"            
## [51] "r_FISHID[42,Intercept]"             "r_FISHID[43,Intercept]"            
## [53] "r_FISHID[44,Intercept]"             "r_FISHID[45,Intercept]"            
## [55] "r_FISHID[46,Intercept]"             "r_FISHID[47,Intercept]"            
## [57] "r_FISHID[48,Intercept]"             "r_FISHID[49,Intercept]"            
## [59] "r_FISHID[50,Intercept]"             "r_FISHID[51,Intercept]"            
## [61] "r_FISHID[52,Intercept]"             "r_FISHID[53,Intercept]"            
## [63] "r_FISHID[54,Intercept]"             "r_FISHID[55,Intercept]"            
## [65] "r_FISHID[56,Intercept]"             "r_FISHID[57,Intercept]"            
## [67] "r_FISHID[58,Intercept]"             "r_FISHID[59,Intercept]"            
## [69] "r_FISHID[60,Intercept]"             "prior_Intercept"                   
## [71] "prior_b_TRIALHypoxia"               "prior_b_TRIALLowSalinity"          
## [73] "prior_b_SMR_contr"                  "prior_b_MASS"                      
## [75] "prior_b_TRIALHypoxia:SMR_contr"     "prior_b_TRIALLowSalinity:SMR_contr"
## [77] "prior_sigma"                        "prior_sd_FISHID"                   
## [79] "lprior"                             "lp__"                              
## [81] "accept_stat__"                      "treedepth__"                       
## [83] "stepsize__"                         "divergent__"                       
## [85] "n_leapfrog__"                       "energy__"

norin.brm3 %>% hypothesis('TRIALHypoxia=0') %>% plot

norin.brm3 %>% hypothesis('SMR_contr=0') %>% plot

## norin.brm3 %>% SUYR_prior_and_posterior()

Random intercept/slope model

While we are here, we might like to explore a random intercept/slope model

priors <- prior(normal(16, 35), class = 'Intercept') +
    prior(normal(0, 70), class = 'b', coef = 'TRIALHypoxia') +
    prior(normal(0, 70), class = 'b', coef = 'TRIALLowSalinity') +
    prior(normal(0, 54), class = 'b', coef = 'SMR_contr') +
    prior(normal(0, 4), class = 'b', coef = 'MASS') +
    ## prior(normal(0, 70), class = 'b') +
    prior(student_t(3,0,34), class = 'sigma') +
    prior(cauchy(0, 5.8), class = 'sd') +
    prior(lkj_corr_cholesky(1), class = 'cor')
norin.form <- bf(CHANGE ~ (TRIAL|FISHID) + TRIAL*SMR_contr + MASS,
                 ## sigma ~ TRIAL*SMR_contr + offset(MASS) + (1|FISHID),
                 family = gaussian()
                 )
norin.brm4 <-  brm(norin.form, 
                  data = norin,
                  prior = priors,
                  sample_prior = 'yes',
                  iter = 5000,
                  warmup = 1000,
                  chains = 3, cores = 3,
                  thin = 5,
                  refresh = 0,
                  control = list(adapt_delta=0.99),
                  backend = "cmdstanr"
                  )

## Start sampling

## Running MCMC with 3 parallel chains...
## 
## Chain 3 finished in 31.1 seconds.
## Chain 1 finished in 42.2 seconds.
## Chain 2 finished in 53.3 seconds.
## 
## All 3 chains finished successfully.
## Mean chain execution time: 42.2 seconds.
## Total execution time: 53.4 seconds.

## Warning: 4 of 2400 (0.0%) transitions hit the maximum treedepth limit of 10.
## See https://mc-stan.org/misc/warnings for details.

save(norin.brm4, file = '../ws/testing/norin.brm4')

norin.brm4 %>% get_variables()

##   [1] "b_Intercept"                               
##   [2] "b_TRIALHypoxia"                            
##   [3] "b_TRIALLowSalinity"                        
##   [4] "b_SMR_contr"                               
##   [5] "b_MASS"                                    
##   [6] "b_TRIALHypoxia:SMR_contr"                  
##   [7] "b_TRIALLowSalinity:SMR_contr"              
##   [8] "sd_FISHID__Intercept"                      
##   [9] "sd_FISHID__TRIALHypoxia"                   
##  [10] "sd_FISHID__TRIALLowSalinity"               
##  [11] "cor_FISHID__Intercept__TRIALHypoxia"       
##  [12] "cor_FISHID__Intercept__TRIALLowSalinity"   
##  [13] "cor_FISHID__TRIALHypoxia__TRIALLowSalinity"
##  [14] "sigma"                                     
##  [15] "r_FISHID[1,Intercept]"                     
##  [16] "r_FISHID[2,Intercept]"                     
##  [17] "r_FISHID[3,Intercept]"                     
##  [18] "r_FISHID[4,Intercept]"                     
##  [19] "r_FISHID[5,Intercept]"                     
##  [20] "r_FISHID[6,Intercept]"                     
##  [21] "r_FISHID[7,Intercept]"                     
##  [22] "r_FISHID[8,Intercept]"                     
##  [23] "r_FISHID[9,Intercept]"                     
##  [24] "r_FISHID[10,Intercept]"                    
##  [25] "r_FISHID[11,Intercept]"                    
##  [26] "r_FISHID[12,Intercept]"                    
##  [27] "r_FISHID[13,Intercept]"                    
##  [28] "r_FISHID[14,Intercept]"                    
##  [29] "r_FISHID[15,Intercept]"                    
##  [30] "r_FISHID[16,Intercept]"                    
##  [31] "r_FISHID[17,Intercept]"                    
##  [32] "r_FISHID[18,Intercept]"                    
##  [33] "r_FISHID[19,Intercept]"                    
##  [34] "r_FISHID[20,Intercept]"                    
##  [35] "r_FISHID[21,Intercept]"                    
##  [36] "r_FISHID[22,Intercept]"                    
##  [37] "r_FISHID[23,Intercept]"                    
##  [38] "r_FISHID[24,Intercept]"                    
##  [39] "r_FISHID[25,Intercept]"                    
##  [40] "r_FISHID[26,Intercept]"                    
##  [41] "r_FISHID[27,Intercept]"                    
##  [42] "r_FISHID[28,Intercept]"                    
##  [43] "r_FISHID[29,Intercept]"                    
##  [44] "r_FISHID[30,Intercept]"                    
##  [45] "r_FISHID[31,Intercept]"                    
##  [46] "r_FISHID[32,Intercept]"                    
##  [47] "r_FISHID[33,Intercept]"                    
##  [48] "r_FISHID[34,Intercept]"                    
##  [49] "r_FISHID[35,Intercept]"                    
##  [50] "r_FISHID[36,Intercept]"                    
##  [51] "r_FISHID[37,Intercept]"                    
##  [52] "r_FISHID[38,Intercept]"                    
##  [53] "r_FISHID[39,Intercept]"                    
##  [54] "r_FISHID[40,Intercept]"                    
##  [55] "r_FISHID[41,Intercept]"                    
##  [56] "r_FISHID[42,Intercept]"                    
##  [57] "r_FISHID[43,Intercept]"                    
##  [58] "r_FISHID[44,Intercept]"                    
##  [59] "r_FISHID[45,Intercept]"                    
##  [60] "r_FISHID[46,Intercept]"                    
##  [61] "r_FISHID[47,Intercept]"                    
##  [62] "r_FISHID[48,Intercept]"                    
##  [63] "r_FISHID[49,Intercept]"                    
##  [64] "r_FISHID[50,Intercept]"                    
##  [65] "r_FISHID[51,Intercept]"                    
##  [66] "r_FISHID[52,Intercept]"                    
##  [67] "r_FISHID[53,Intercept]"                    
##  [68] "r_FISHID[54,Intercept]"                    
##  [69] "r_FISHID[55,Intercept]"                    
##  [70] "r_FISHID[56,Intercept]"                    
##  [71] "r_FISHID[57,Intercept]"                    
##  [72] "r_FISHID[58,Intercept]"                    
##  [73] "r_FISHID[59,Intercept]"                    
##  [74] "r_FISHID[60,Intercept]"                    
##  [75] "r_FISHID[1,TRIALHypoxia]"                  
##  [76] "r_FISHID[2,TRIALHypoxia]"                  
##  [77] "r_FISHID[3,TRIALHypoxia]"                  
##  [78] "r_FISHID[4,TRIALHypoxia]"                  
##  [79] "r_FISHID[5,TRIALHypoxia]"                  
##  [80] "r_FISHID[6,TRIALHypoxia]"                  
##  [81] "r_FISHID[7,TRIALHypoxia]"                  
##  [82] "r_FISHID[8,TRIALHypoxia]"                  
##  [83] "r_FISHID[9,TRIALHypoxia]"                  
##  [84] "r_FISHID[10,TRIALHypoxia]"                 
##  [85] "r_FISHID[11,TRIALHypoxia]"                 
##  [86] "r_FISHID[12,TRIALHypoxia]"                 
##  [87] "r_FISHID[13,TRIALHypoxia]"                 
##  [88] "r_FISHID[14,TRIALHypoxia]"                 
##  [89] "r_FISHID[15,TRIALHypoxia]"                 
##  [90] "r_FISHID[16,TRIALHypoxia]"                 
##  [91] "r_FISHID[17,TRIALHypoxia]"                 
##  [92] "r_FISHID[18,TRIALHypoxia]"                 
##  [93] "r_FISHID[19,TRIALHypoxia]"                 
##  [94] "r_FISHID[20,TRIALHypoxia]"                 
##  [95] "r_FISHID[21,TRIALHypoxia]"                 
##  [96] "r_FISHID[22,TRIALHypoxia]"                 
##  [97] "r_FISHID[23,TRIALHypoxia]"                 
##  [98] "r_FISHID[24,TRIALHypoxia]"                 
##  [99] "r_FISHID[25,TRIALHypoxia]"                 
## [100] "r_FISHID[26,TRIALHypoxia]"                 
## [101] "r_FISHID[27,TRIALHypoxia]"                 
## [102] "r_FISHID[28,TRIALHypoxia]"                 
## [103] "r_FISHID[29,TRIALHypoxia]"                 
## [104] "r_FISHID[30,TRIALHypoxia]"                 
## [105] "r_FISHID[31,TRIALHypoxia]"                 
## [106] "r_FISHID[32,TRIALHypoxia]"                 
## [107] "r_FISHID[33,TRIALHypoxia]"                 
## [108] "r_FISHID[34,TRIALHypoxia]"                 
## [109] "r_FISHID[35,TRIALHypoxia]"                 
## [110] "r_FISHID[36,TRIALHypoxia]"                 
## [111] "r_FISHID[37,TRIALHypoxia]"                 
## [112] "r_FISHID[38,TRIALHypoxia]"                 
## [113] "r_FISHID[39,TRIALHypoxia]"                 
## [114] "r_FISHID[40,TRIALHypoxia]"                 
## [115] "r_FISHID[41,TRIALHypoxia]"                 
## [116] "r_FISHID[42,TRIALHypoxia]"                 
## [117] "r_FISHID[43,TRIALHypoxia]"                 
## [118] "r_FISHID[44,TRIALHypoxia]"                 
## [119] "r_FISHID[45,TRIALHypoxia]"                 
## [120] "r_FISHID[46,TRIALHypoxia]"                 
## [121] "r_FISHID[47,TRIALHypoxia]"                 
## [122] "r_FISHID[48,TRIALHypoxia]"                 
## [123] "r_FISHID[49,TRIALHypoxia]"                 
## [124] "r_FISHID[50,TRIALHypoxia]"                 
## [125] "r_FISHID[51,TRIALHypoxia]"                 
## [126] "r_FISHID[52,TRIALHypoxia]"                 
## [127] "r_FISHID[53,TRIALHypoxia]"                 
## [128] "r_FISHID[54,TRIALHypoxia]"                 
## [129] "r_FISHID[55,TRIALHypoxia]"                 
## [130] "r_FISHID[56,TRIALHypoxia]"                 
## [131] "r_FISHID[57,TRIALHypoxia]"                 
## [132] "r_FISHID[58,TRIALHypoxia]"                 
## [133] "r_FISHID[59,TRIALHypoxia]"                 
## [134] "r_FISHID[60,TRIALHypoxia]"                 
## [135] "r_FISHID[1,TRIALLowSalinity]"              
## [136] "r_FISHID[2,TRIALLowSalinity]"              
## [137] "r_FISHID[3,TRIALLowSalinity]"              
## [138] "r_FISHID[4,TRIALLowSalinity]"              
## [139] "r_FISHID[5,TRIALLowSalinity]"              
## [140] "r_FISHID[6,TRIALLowSalinity]"              
## [141] "r_FISHID[7,TRIALLowSalinity]"              
## [142] "r_FISHID[8,TRIALLowSalinity]"              
## [143] "r_FISHID[9,TRIALLowSalinity]"              
## [144] "r_FISHID[10,TRIALLowSalinity]"             
## [145] "r_FISHID[11,TRIALLowSalinity]"             
## [146] "r_FISHID[12,TRIALLowSalinity]"             
## [147] "r_FISHID[13,TRIALLowSalinity]"             
## [148] "r_FISHID[14,TRIALLowSalinity]"             
## [149] "r_FISHID[15,TRIALLowSalinity]"             
## [150] "r_FISHID[16,TRIALLowSalinity]"             
## [151] "r_FISHID[17,TRIALLowSalinity]"             
## [152] "r_FISHID[18,TRIALLowSalinity]"             
## [153] "r_FISHID[19,TRIALLowSalinity]"             
## [154] "r_FISHID[20,TRIALLowSalinity]"             
## [155] "r_FISHID[21,TRIALLowSalinity]"             
## [156] "r_FISHID[22,TRIALLowSalinity]"             
## [157] "r_FISHID[23,TRIALLowSalinity]"             
## [158] "r_FISHID[24,TRIALLowSalinity]"             
## [159] "r_FISHID[25,TRIALLowSalinity]"             
## [160] "r_FISHID[26,TRIALLowSalinity]"             
## [161] "r_FISHID[27,TRIALLowSalinity]"             
## [162] "r_FISHID[28,TRIALLowSalinity]"             
## [163] "r_FISHID[29,TRIALLowSalinity]"             
## [164] "r_FISHID[30,TRIALLowSalinity]"             
## [165] "r_FISHID[31,TRIALLowSalinity]"             
## [166] "r_FISHID[32,TRIALLowSalinity]"             
## [167] "r_FISHID[33,TRIALLowSalinity]"             
## [168] "r_FISHID[34,TRIALLowSalinity]"             
## [169] "r_FISHID[35,TRIALLowSalinity]"             
## [170] "r_FISHID[36,TRIALLowSalinity]"             
## [171] "r_FISHID[37,TRIALLowSalinity]"             
## [172] "r_FISHID[38,TRIALLowSalinity]"             
## [173] "r_FISHID[39,TRIALLowSalinity]"             
## [174] "r_FISHID[40,TRIALLowSalinity]"             
## [175] "r_FISHID[41,TRIALLowSalinity]"             
## [176] "r_FISHID[42,TRIALLowSalinity]"             
## [177] "r_FISHID[43,TRIALLowSalinity]"             
## [178] "r_FISHID[44,TRIALLowSalinity]"             
## [179] "r_FISHID[45,TRIALLowSalinity]"             
## [180] "r_FISHID[46,TRIALLowSalinity]"             
## [181] "r_FISHID[47,TRIALLowSalinity]"             
## [182] "r_FISHID[48,TRIALLowSalinity]"             
## [183] "r_FISHID[49,TRIALLowSalinity]"             
## [184] "r_FISHID[50,TRIALLowSalinity]"             
## [185] "r_FISHID[51,TRIALLowSalinity]"             
## [186] "r_FISHID[52,TRIALLowSalinity]"             
## [187] "r_FISHID[53,TRIALLowSalinity]"             
## [188] "r_FISHID[54,TRIALLowSalinity]"             
## [189] "r_FISHID[55,TRIALLowSalinity]"             
## [190] "r_FISHID[56,TRIALLowSalinity]"             
## [191] "r_FISHID[57,TRIALLowSalinity]"             
## [192] "r_FISHID[58,TRIALLowSalinity]"             
## [193] "r_FISHID[59,TRIALLowSalinity]"             
## [194] "r_FISHID[60,TRIALLowSalinity]"             
## [195] "prior_Intercept"                           
## [196] "prior_b_TRIALHypoxia"                      
## [197] "prior_b_TRIALLowSalinity"                  
## [198] "prior_b_SMR_contr"                         
## [199] "prior_b_MASS"                              
## [200] "prior_sigma"                               
## [201] "prior_sd_FISHID"                           
## [202] "prior_cor_FISHID"                          
## [203] "lprior"                                    
## [204] "lp__"                                      
## [205] "accept_stat__"                             
## [206] "treedepth__"                               
## [207] "stepsize__"                                
## [208] "divergent__"                               
## [209] "n_leapfrog__"                              
## [210] "energy__"

norin.brm4 %>% hypothesis('TRIALHypoxia=0') %>% plot

norin.brm4 %>% hypothesis('SMR_contr=0') %>% plot

## norin.brm4 %>% SUYR_prior_and_posterior()

norin.brm4 %>%
  posterior_samples %>%
  dplyr::select(-`lp__`) %>%
  pivot_longer(everything(), names_to = 'key') %>% 
  filter(!str_detect(key, '^r')) %>%
  mutate(Type = ifelse(str_detect(key, 'prior'), 'Prior', 'Posterior'),
         Class = case_when(
             str_detect(key, '(^b|^prior).*Intercept$') ~ 'Intercept',
             str_detect(key, 'b_TRIAL.*|prior_b_TRIAL.*') & !str_detect(key, '.*\\:.*') ~ 'TRIAL',
             str_detect(key, 'b_SMR_contr|prior_b_SMR_contr') ~ 'SMR',
             str_detect(key, 'b_MASS|prior_b_MASS') ~ 'MASS',
             str_detect(key, '.*\\:.*|prior_b_.*\\:.*') ~ 'Interaction',
             str_detect(key, 'sd') ~ 'sd',
             str_detect(key, '^cor|prior_cor') ~ 'cor',
             str_detect(key, 'sigma') ~ 'sigma'),
         Par = str_replace(key, 'b_', '')) %>%
  ggplot(aes(x = Type,  y = value, color = Par)) +
  stat_pointinterval(position = position_dodge())+
  facet_wrap(~Class,  scales = 'free')

## Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
## recommended alternatives.

Compare models

We can compare the two models using LOO

(l.1 <- norin.brm3 %>% loo()) 

## 
## Computed from 2400 by 180 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo   -831.8 13.6
## p_loo        23.8  3.8
## looic      1663.7 27.3
## ------
## Monte Carlo SE of elpd_loo is 0.2.
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. n_eff
## (-Inf, 0.5]   (good)     171   95.0%   205       
##  (0.5, 0.7]   (ok)         9    5.0%   191       
##    (0.7, 1]   (bad)        0    0.0%   <NA>      
##    (1, Inf)   (very bad)   0    0.0%   <NA>      
## 
## All Pareto k estimates are ok (k < 0.7).
## See help('pareto-k-diagnostic') for details.

(l.2 <- norin.brm4 %>% loo())

## Warning: Found 56 observations with a pareto_k > 0.7 in model '.'. It is
## recommended to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.

## 
## Computed from 2400 by 180 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo   -782.0  9.4
## p_loo        89.8  6.5
## looic      1564.0 18.8
## ------
## Monte Carlo SE of elpd_loo is NA.
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. n_eff
## (-Inf, 0.5]   (good)     55    30.6%   108       
##  (0.5, 0.7]   (ok)       69    38.3%   20        
##    (0.7, 1]   (bad)      49    27.2%   5         
##    (1, Inf)   (very bad)  7     3.9%   2         
## See help('pareto-k-diagnostic') for details.

loo_compare(l.1, l.2)

##   elpd_diff se_diff
## .   0.0       0.0  
## . -49.8       9.4

There is substantially more support for the more complex random intercept/slope model over the simpler random intercept only model. Consequently, we will continue with the random intercept/slope model.

MCMC sampling diagnostics

In addition to the regular model diagnostics checking, for Bayesian analyses, it is also necessary to explore the MCMC sampling diagnostics to be sure that the chains are well mixed and have converged on a stable posterior.

There are a wide variety of tests that range from the big picture, overal chain characteristics to the very specific detailed tests that allow the experienced modeller to drill down to the very fine details of the chain behaviour. Furthermore, there are a multitude of packages and approaches for exploring these diagnostics.

brms

bayesplot

The bayesplot package offers a range of MCMC diagnostics as well as Posterior Probability Checks (PPC), all of which have a convenient plot() interface. Lets start with the MCMC diagnostics.

See list of available diagnostics by name

available_mcmc()

## bayesplot MCMC module:
##   mcmc_acf
##   mcmc_acf_bar
##   mcmc_areas
##   mcmc_areas_ridges
##   mcmc_combo
##   mcmc_dens
##   mcmc_dens_chains
##   mcmc_dens_overlay
##   mcmc_hex
##   mcmc_hist
##   mcmc_hist_by_chain
##   mcmc_intervals
##   mcmc_neff
##   mcmc_neff_hist
##   mcmc_nuts_acceptance
##   mcmc_nuts_divergence
##   mcmc_nuts_energy
##   mcmc_nuts_stepsize
##   mcmc_nuts_treedepth
##   mcmc_pairs
##   mcmc_parcoord
##   mcmc_rank_hist
##   mcmc_rank_overlay
##   mcmc_recover_hist
##   mcmc_recover_intervals
##   mcmc_recover_scatter
##   mcmc_rhat
##   mcmc_rhat_hist
##   mcmc_scatter
##   mcmc_trace
##   mcmc_trace_highlight
##   mcmc_violin

Of these, we will focus on:

• trace: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different shade of blue, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).

norin.brm4 %>% mcmc_plot(type='trace')

## No divergences to plot.

The chains appear well mixed and very similar

• acf_bar (autocorrelation function): plots the autocorrelation between successive MCMC sample lags for each parameter and each chain

norin.brm4 %>% mcmc_plot(type='acf_bar')

There IS evidence of autocorrelation in the MCMC samples. We should consider thinning further!

• rhat_hist: Rhat is a scale reduction factor measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.

norin.brm4 %>% mcmc_plot(type='rhat_hist')

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

All Rhat values are below 1.05, suggesting the chains have converged.

• neff_hist (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).

  If the ratios are low, tightening the priors may help.

norin.brm4 %>% mcmc_plot(type='neff_hist')

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

There are a small number of parameters that have very low effective sample sizes. It might be worth exploring the cause(s) of this to determine whether it is concerning.

More diagnostics

norin.brm4 %>% mcmc_plot(type='combo')

norin.brm4 %>% mcmc_plot(type='violin')

stan plots

The rstan package offers a range of MCMC diagnostics. Lets start with the MCMC diagnostics.

Of these, we will focus on:

• stan_trace: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different colour, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).

norin.brm4 %>% get_variables()

##   [1] "b_Intercept"                               
##   [2] "b_TRIALHypoxia"                            
##   [3] "b_TRIALLowSalinity"                        
##   [4] "b_SMR_contr"                               
##   [5] "b_MASS"                                    
##   [6] "b_TRIALHypoxia:SMR_contr"                  
##   [7] "b_TRIALLowSalinity:SMR_contr"              
##   [8] "sd_FISHID__Intercept"                      
##   [9] "sd_FISHID__TRIALHypoxia"                   
##  [10] "sd_FISHID__TRIALLowSalinity"               
##  [11] "cor_FISHID__Intercept__TRIALHypoxia"       
##  [12] "cor_FISHID__Intercept__TRIALLowSalinity"   
##  [13] "cor_FISHID__TRIALHypoxia__TRIALLowSalinity"
##  [14] "sigma"                                     
##  [15] "r_FISHID[1,Intercept]"                     
##  [16] "r_FISHID[2,Intercept]"                     
##  [17] "r_FISHID[3,Intercept]"                     
##  [18] "r_FISHID[4,Intercept]"                     
##  [19] "r_FISHID[5,Intercept]"                     
##  [20] "r_FISHID[6,Intercept]"                     
##  [21] "r_FISHID[7,Intercept]"                     
##  [22] "r_FISHID[8,Intercept]"                     
##  [23] "r_FISHID[9,Intercept]"                     
##  [24] "r_FISHID[10,Intercept]"                    
##  [25] "r_FISHID[11,Intercept]"                    
##  [26] "r_FISHID[12,Intercept]"                    
##  [27] "r_FISHID[13,Intercept]"                    
##  [28] "r_FISHID[14,Intercept]"                    
##  [29] "r_FISHID[15,Intercept]"                    
##  [30] "r_FISHID[16,Intercept]"                    
##  [31] "r_FISHID[17,Intercept]"                    
##  [32] "r_FISHID[18,Intercept]"                    
##  [33] "r_FISHID[19,Intercept]"                    
##  [34] "r_FISHID[20,Intercept]"                    
##  [35] "r_FISHID[21,Intercept]"                    
##  [36] "r_FISHID[22,Intercept]"                    
##  [37] "r_FISHID[23,Intercept]"                    
##  [38] "r_FISHID[24,Intercept]"                    
##  [39] "r_FISHID[25,Intercept]"                    
##  [40] "r_FISHID[26,Intercept]"                    
##  [41] "r_FISHID[27,Intercept]"                    
##  [42] "r_FISHID[28,Intercept]"                    
##  [43] "r_FISHID[29,Intercept]"                    
##  [44] "r_FISHID[30,Intercept]"                    
##  [45] "r_FISHID[31,Intercept]"                    
##  [46] "r_FISHID[32,Intercept]"                    
##  [47] "r_FISHID[33,Intercept]"                    
##  [48] "r_FISHID[34,Intercept]"                    
##  [49] "r_FISHID[35,Intercept]"                    
##  [50] "r_FISHID[36,Intercept]"                    
##  [51] "r_FISHID[37,Intercept]"                    
##  [52] "r_FISHID[38,Intercept]"                    
##  [53] "r_FISHID[39,Intercept]"                    
##  [54] "r_FISHID[40,Intercept]"                    
##  [55] "r_FISHID[41,Intercept]"                    
##  [56] "r_FISHID[42,Intercept]"                    
##  [57] "r_FISHID[43,Intercept]"                    
##  [58] "r_FISHID[44,Intercept]"                    
##  [59] "r_FISHID[45,Intercept]"                    
##  [60] "r_FISHID[46,Intercept]"                    
##  [61] "r_FISHID[47,Intercept]"                    
##  [62] "r_FISHID[48,Intercept]"                    
##  [63] "r_FISHID[49,Intercept]"                    
##  [64] "r_FISHID[50,Intercept]"                    
##  [65] "r_FISHID[51,Intercept]"                    
##  [66] "r_FISHID[52,Intercept]"                    
##  [67] "r_FISHID[53,Intercept]"                    
##  [68] "r_FISHID[54,Intercept]"                    
##  [69] "r_FISHID[55,Intercept]"                    
##  [70] "r_FISHID[56,Intercept]"                    
##  [71] "r_FISHID[57,Intercept]"                    
##  [72] "r_FISHID[58,Intercept]"                    
##  [73] "r_FISHID[59,Intercept]"                    
##  [74] "r_FISHID[60,Intercept]"                    
##  [75] "r_FISHID[1,TRIALHypoxia]"                  
##  [76] "r_FISHID[2,TRIALHypoxia]"                  
##  [77] "r_FISHID[3,TRIALHypoxia]"                  
##  [78] "r_FISHID[4,TRIALHypoxia]"                  
##  [79] "r_FISHID[5,TRIALHypoxia]"                  
##  [80] "r_FISHID[6,TRIALHypoxia]"                  
##  [81] "r_FISHID[7,TRIALHypoxia]"                  
##  [82] "r_FISHID[8,TRIALHypoxia]"                  
##  [83] "r_FISHID[9,TRIALHypoxia]"                  
##  [84] "r_FISHID[10,TRIALHypoxia]"                 
##  [85] "r_FISHID[11,TRIALHypoxia]"                 
##  [86] "r_FISHID[12,TRIALHypoxia]"                 
##  [87] "r_FISHID[13,TRIALHypoxia]"                 
##  [88] "r_FISHID[14,TRIALHypoxia]"                 
##  [89] "r_FISHID[15,TRIALHypoxia]"                 
##  [90] "r_FISHID[16,TRIALHypoxia]"                 
##  [91] "r_FISHID[17,TRIALHypoxia]"                 
##  [92] "r_FISHID[18,TRIALHypoxia]"                 
##  [93] "r_FISHID[19,TRIALHypoxia]"                 
##  [94] "r_FISHID[20,TRIALHypoxia]"                 
##  [95] "r_FISHID[21,TRIALHypoxia]"                 
##  [96] "r_FISHID[22,TRIALHypoxia]"                 
##  [97] "r_FISHID[23,TRIALHypoxia]"                 
##  [98] "r_FISHID[24,TRIALHypoxia]"                 
##  [99] "r_FISHID[25,TRIALHypoxia]"                 
## [100] "r_FISHID[26,TRIALHypoxia]"                 
## [101] "r_FISHID[27,TRIALHypoxia]"                 
## [102] "r_FISHID[28,TRIALHypoxia]"                 
## [103] "r_FISHID[29,TRIALHypoxia]"                 
## [104] "r_FISHID[30,TRIALHypoxia]"                 
## [105] "r_FISHID[31,TRIALHypoxia]"                 
## [106] "r_FISHID[32,TRIALHypoxia]"                 
## [107] "r_FISHID[33,TRIALHypoxia]"                 
## [108] "r_FISHID[34,TRIALHypoxia]"                 
## [109] "r_FISHID[35,TRIALHypoxia]"                 
## [110] "r_FISHID[36,TRIALHypoxia]"                 
## [111] "r_FISHID[37,TRIALHypoxia]"                 
## [112] "r_FISHID[38,TRIALHypoxia]"                 
## [113] "r_FISHID[39,TRIALHypoxia]"                 
## [114] "r_FISHID[40,TRIALHypoxia]"                 
## [115] "r_FISHID[41,TRIALHypoxia]"                 
## [116] "r_FISHID[42,TRIALHypoxia]"                 
## [117] "r_FISHID[43,TRIALHypoxia]"                 
## [118] "r_FISHID[44,TRIALHypoxia]"                 
## [119] "r_FISHID[45,TRIALHypoxia]"                 
## [120] "r_FISHID[46,TRIALHypoxia]"                 
## [121] "r_FISHID[47,TRIALHypoxia]"                 
## [122] "r_FISHID[48,TRIALHypoxia]"                 
## [123] "r_FISHID[49,TRIALHypoxia]"                 
## [124] "r_FISHID[50,TRIALHypoxia]"                 
## [125] "r_FISHID[51,TRIALHypoxia]"                 
## [126] "r_FISHID[52,TRIALHypoxia]"                 
## [127] "r_FISHID[53,TRIALHypoxia]"                 
## [128] "r_FISHID[54,TRIALHypoxia]"                 
## [129] "r_FISHID[55,TRIALHypoxia]"                 
## [130] "r_FISHID[56,TRIALHypoxia]"                 
## [131] "r_FISHID[57,TRIALHypoxia]"                 
## [132] "r_FISHID[58,TRIALHypoxia]"                 
## [133] "r_FISHID[59,TRIALHypoxia]"                 
## [134] "r_FISHID[60,TRIALHypoxia]"                 
## [135] "r_FISHID[1,TRIALLowSalinity]"              
## [136] "r_FISHID[2,TRIALLowSalinity]"              
## [137] "r_FISHID[3,TRIALLowSalinity]"              
## [138] "r_FISHID[4,TRIALLowSalinity]"              
## [139] "r_FISHID[5,TRIALLowSalinity]"              
## [140] "r_FISHID[6,TRIALLowSalinity]"              
## [141] "r_FISHID[7,TRIALLowSalinity]"              
## [142] "r_FISHID[8,TRIALLowSalinity]"              
## [143] "r_FISHID[9,TRIALLowSalinity]"              
## [144] "r_FISHID[10,TRIALLowSalinity]"             
## [145] "r_FISHID[11,TRIALLowSalinity]"             
## [146] "r_FISHID[12,TRIALLowSalinity]"             
## [147] "r_FISHID[13,TRIALLowSalinity]"             
## [148] "r_FISHID[14,TRIALLowSalinity]"             
## [149] "r_FISHID[15,TRIALLowSalinity]"             
## [150] "r_FISHID[16,TRIALLowSalinity]"             
## [151] "r_FISHID[17,TRIALLowSalinity]"             
## [152] "r_FISHID[18,TRIALLowSalinity]"             
## [153] "r_FISHID[19,TRIALLowSalinity]"             
## [154] "r_FISHID[20,TRIALLowSalinity]"             
## [155] "r_FISHID[21,TRIALLowSalinity]"             
## [156] "r_FISHID[22,TRIALLowSalinity]"             
## [157] "r_FISHID[23,TRIALLowSalinity]"             
## [158] "r_FISHID[24,TRIALLowSalinity]"             
## [159] "r_FISHID[25,TRIALLowSalinity]"             
## [160] "r_FISHID[26,TRIALLowSalinity]"             
## [161] "r_FISHID[27,TRIALLowSalinity]"             
## [162] "r_FISHID[28,TRIALLowSalinity]"             
## [163] "r_FISHID[29,TRIALLowSalinity]"             
## [164] "r_FISHID[30,TRIALLowSalinity]"             
## [165] "r_FISHID[31,TRIALLowSalinity]"             
## [166] "r_FISHID[32,TRIALLowSalinity]"             
## [167] "r_FISHID[33,TRIALLowSalinity]"             
## [168] "r_FISHID[34,TRIALLowSalinity]"             
## [169] "r_FISHID[35,TRIALLowSalinity]"             
## [170] "r_FISHID[36,TRIALLowSalinity]"             
## [171] "r_FISHID[37,TRIALLowSalinity]"             
## [172] "r_FISHID[38,TRIALLowSalinity]"             
## [173] "r_FISHID[39,TRIALLowSalinity]"             
## [174] "r_FISHID[40,TRIALLowSalinity]"             
## [175] "r_FISHID[41,TRIALLowSalinity]"             
## [176] "r_FISHID[42,TRIALLowSalinity]"             
## [177] "r_FISHID[43,TRIALLowSalinity]"             
## [178] "r_FISHID[44,TRIALLowSalinity]"             
## [179] "r_FISHID[45,TRIALLowSalinity]"             
## [180] "r_FISHID[46,TRIALLowSalinity]"             
## [181] "r_FISHID[47,TRIALLowSalinity]"             
## [182] "r_FISHID[48,TRIALLowSalinity]"             
## [183] "r_FISHID[49,TRIALLowSalinity]"             
## [184] "r_FISHID[50,TRIALLowSalinity]"             
## [185] "r_FISHID[51,TRIALLowSalinity]"             
## [186] "r_FISHID[52,TRIALLowSalinity]"             
## [187] "r_FISHID[53,TRIALLowSalinity]"             
## [188] "r_FISHID[54,TRIALLowSalinity]"             
## [189] "r_FISHID[55,TRIALLowSalinity]"             
## [190] "r_FISHID[56,TRIALLowSalinity]"             
## [191] "r_FISHID[57,TRIALLowSalinity]"             
## [192] "r_FISHID[58,TRIALLowSalinity]"             
## [193] "r_FISHID[59,TRIALLowSalinity]"             
## [194] "r_FISHID[60,TRIALLowSalinity]"             
## [195] "prior_Intercept"                           
## [196] "prior_b_TRIALHypoxia"                      
## [197] "prior_b_TRIALLowSalinity"                  
## [198] "prior_b_SMR_contr"                         
## [199] "prior_b_MASS"                              
## [200] "prior_sigma"                               
## [201] "prior_sd_FISHID"                           
## [202] "prior_cor_FISHID"                          
## [203] "lprior"                                    
## [204] "lp__"                                      
## [205] "accept_stat__"                             
## [206] "treedepth__"                               
## [207] "stepsize__"                                
## [208] "divergent__"                               
## [209] "n_leapfrog__"                              
## [210] "energy__"

pars <- norin.brm4 %>% get_variables()
pars <- str_extract(pars, '^b_.*|^sigma$|^sd.*') %>% na.omit()

norin.brm4$fit %>%
    stan_trace(pars = pars)

The chains appear well mixed and very similar

• stan_acf (autocorrelation function): plots the autocorrelation between successive MCMC sample lags for each parameter and each chain

norin.brm4$fit %>%
    stan_ac(pars = pars)

There IS evidence of autocorrelation in the MCMC samples. We should consider thinning further!

• stan_rhat: Rhat is a scale reduction factor measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.

norin.brm4$fit %>% stan_rhat() 

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

All Rhat values are below 1.05, suggesting the chains have converged.

• stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).

  If the ratios are low, tightening the priors may help.

norin.brm4$fit %>% stan_ess()

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

There are a small number of parameters that have very low effective sample sizes. It might be worth exploring the cause(s) of this to determine whether it is concerning.

norin.brm4$fit %>%
    stan_dens(separate_chains = TRUE, pars = pars)

ggmcmc

The ggmean package also as a set of MCMC diagnostic functions. Lets start with the MCMC diagnostics.

Of these, we will focus on:

• ggs_traceplot: this plots the estimates of each parameter over the post-warmup length of each MCMC chain. Each chain is plotted in a different colour, with each parameter in its own facet. Ideally, each trace should just look like noise without any discernible drift and each of the traces for a specific parameter should look the same (i.e, should not be displaced above or below any other trace for that parameter).

norin.ggs <- norin.brm4 %>% ggs(burnin = FALSE, inc_warmup = FALSE)

## Warning in custom.sort(D$Parameter): NAs introduced by coercion

norin.ggs %>% ggs_traceplot()

The chains appear well mixed and very similar

• gss_autocorrelation (autocorrelation function): plots the autocorrelation between successive MCMC sample lags for each parameter and each chain

ggs_autocorrelation(norin.ggs)

There is no evidence of autocorrelation in the MCMC samples

• stan_rhat: Rhat is a scale reduction factor measure of convergence between the chains. The closer the values are to 1, the more the chains have converged. Values greater than 1.05 indicate a lack of convergence. There will be an Rhat value for each parameter estimated.

ggs_Rhat(norin.ggs, scaling = 1.01)

All Rhat values are below 1.01, suggesting the chains have converged.

• stan_ess (number of effective samples): the ratio of the number of effective samples (those not rejected by the sampler) to the number of samples provides an indication of the effectiveness (and efficiency) of the MCMC sampler. Ratios that are less than 0.5 for a parameter suggest that the sampler spent considerable time in difficult areas of the sampling domain and rejected more than half of the samples (replacing them with the previous effective sample).

  If the ratios are low, tightening the priors may help.

ggs_effective(norin.ggs)

Ratios all very high.

More diagnostics

ggs_crosscorrelation(norin.ggs)

ggs_grb(norin.ggs)

pp check

Post predictive checks provide additional diagnostics about the fit of the model. Specifically, they provide a comparison between predictions drawn from the model and the observed data used to train the model.

See list of available diagnostics by name

available_ppc()

## bayesplot PPC module:
##   ppc_bars
##   ppc_bars_grouped
##   ppc_boxplot
##   ppc_dens
##   ppc_dens_overlay
##   ppc_dens_overlay_grouped
##   ppc_ecdf_overlay
##   ppc_ecdf_overlay_grouped
##   ppc_error_binned
##   ppc_error_hist
##   ppc_error_hist_grouped
##   ppc_error_scatter
##   ppc_error_scatter_avg
##   ppc_error_scatter_avg_grouped
##   ppc_error_scatter_avg_vs_x
##   ppc_freqpoly
##   ppc_freqpoly_grouped
##   ppc_hist
##   ppc_intervals
##   ppc_intervals_grouped
##   ppc_km_overlay
##   ppc_km_overlay_grouped
##   ppc_loo_intervals
##   ppc_loo_pit
##   ppc_loo_pit_overlay
##   ppc_loo_pit_qq
##   ppc_loo_ribbon
##   ppc_ribbon
##   ppc_ribbon_grouped
##   ppc_rootogram
##   ppc_scatter
##   ppc_scatter_avg
##   ppc_scatter_avg_grouped
##   ppc_stat
##   ppc_stat_2d
##   ppc_stat_freqpoly
##   ppc_stat_freqpoly_grouped
##   ppc_stat_grouped
##   ppc_violin_grouped

• dens_overlay: plots the density distribution of the observed data (black line) overlayed ontop of 50 density distributions generated from draws from the model (light blue). Ideally, the 50 realisations should be roughly consistent with the observed data.

norin.brm4 %>% pp_check(type = 'dens_overlay', ndraws = 100)

The model draws appear to be consistent with the observed data.

• error_scatter_avg: this plots the observed values against the average residuals. Similar to a residual plot, we do not want to see any patterns in this plot. Note, this is not really that useful for models that involve a binomial response

norin.brm4 %>% pp_check(type = 'error_scatter_avg')

## Using all posterior draws for ppc type 'error_scatter_avg' by default.

This is not really interpretable

• intervals: plots the observed data overlayed ontop of posterior predictions associated with each level of the predictor. Ideally, the observed data should all fall within the predictive intervals.

norin.brm4 %>% pp_check(group = 'TRIAL', type = 'intervals')

## Using all posterior draws for ppc type 'intervals' by default.

## Warning: The following arguments were unrecognized and ignored: group

norin.brm3 %>% pp_check(group = 'TRIAL', type = 'intervals_grouped')

## Using all posterior draws for ppc type 'intervals_grouped' by default.

norin.brm3 %>% pp_check(group = 'TRIAL', type = 'violin_grouped')

## Using all posterior draws for ppc type 'violin_grouped' by default.

The shinystan package allows the full suite of MCMC diagnostics and posterior predictive checks to be accessed via a web interface.

#library(shinystan)
#launch_shinystan(norin.brm2)

DHARMa residuals

DHARMa residuals provide very useful diagnostics. Unfortunately, we cannot directly use the simulateResiduals() function to generate the simulated residuals. However, if we are willing to calculate some of the components ourself, we can still obtain the simulated residuals from the fitted stan model.

We need to supply:

• simulated (predicted) responses associated with each observation.
• observed values
• fitted (predicted) responses (averaged) associated with each observation

preds <- norin.brm4 %>% posterior_predict(ndraws = 250,  summary = FALSE)
norin.resids <- createDHARMa(simulatedResponse = t(preds),
                            observedResponse = norin$CHANGE,
                            fittedPredictedResponse = apply(preds, 2, median),
                            integerResponse = FALSE)
plot(norin.resids, quantreg = FALSE)

norin.resids %>% testDispersion()

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.31164, p-value < 2.2e-16
## alternative hypothesis: two.sided

Conclusions:

• the simulated residuals do suggest some issues with the residuals
• the model appears to be under-dispersed
• there is evidence of a lack of fit.

Partial effects plots

brms

conditional_effects

norin.brm4 %>%
    conditional_effects() %>%
    plot(ask = FALSE, points = TRUE, plot = FALSE) %>%
    wrap_plots()

ggpredict

norin.brm4 %>%
    ggpredict(terms = c("SMR_contr", "TRIAL")) %>%
    plot(add.data = TRUE)

## Note: uncertainty of error terms are not taken into account. You may want to use `rstantools::posterior_predict()`.

ggemmeans

norin.brm4 %>%
    ggemmeans(~SMR_contr*TRIAL) %>%
    plot(add.data = TRUE)

fitted_draws

It is not really possible to do this via the fitted draws as it would not be marginalising over MASS or FISHID.

Model investigation

brms

brms captures the MCMC samples from stan within the returned list. There are numerous ways to retrieve and summarise these samples. The first three provide convenient numeric summaries from which you can draw conclusions, the last four provide ways of obtaining the full posteriors.

summary

The summary() method generates simple summaries (mean, standard deviation as well as 10, 50 and 90 percentiles).

norin.brm4 %>% summary()

##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: CHANGE ~ (TRIAL | FISHID) + TRIAL * SMR_contr + MASS 
##    Data: norin (Number of observations: 180) 
##   Draws: 3 chains, each with iter = 5000; warmup = 1000; thin = 5;
##          total post-warmup draws = 2400
## 
## Group-Level Effects: 
## ~FISHID (Number of levels: 60) 
##                                    Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept)                         11.85      2.99     5.76    17.61 1.00
## sd(TRIALHypoxia)                       6.59      3.80     0.54    15.13 1.00
## sd(TRIALLowSalinity)                  28.13      3.84    20.90    35.99 1.00
## cor(Intercept,TRIALHypoxia)           -0.32      0.37    -0.86     0.63 1.00
## cor(Intercept,TRIALLowSalinity)        0.07      0.24    -0.33     0.59 1.00
## cor(TRIALHypoxia,TRIALLowSalinity)    -0.51      0.32    -0.95     0.24 1.00
##                                    Bulk_ESS Tail_ESS
## sd(Intercept)                           570      967
## sd(TRIALHypoxia)                        638      457
## sd(TRIALLowSalinity)                   1444     2012
## cor(Intercept,TRIALHypoxia)            2034     1716
## cor(Intercept,TRIALLowSalinity)        1267     1755
## cor(TRIALHypoxia,TRIALLowSalinity)      533      907
## 
## Population-Level Effects: 
##                            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept                    147.67     22.71   103.60   191.43 1.00     2303
## TRIALHypoxia                -107.50     20.87  -148.28   -65.78 1.00     2105
## TRIALLowSalinity             -72.54     32.46  -136.82    -7.76 1.00     2327
## SMR_contr                    -19.70      3.75   -26.98   -12.56 1.00     2332
## MASS                           0.12      0.24    -0.34     0.60 1.00     2547
## TRIALHypoxia:SMR_contr        10.94      4.03     2.82    18.68 1.00     2128
## TRIALLowSalinity:SMR_contr     5.95      6.23    -6.39    17.99 1.00     2354
##                            Tail_ESS
## Intercept                      2345
## TRIALHypoxia                   1997
## TRIALLowSalinity               2370
## SMR_contr                      2200
## MASS                           2334
## TRIALHypoxia:SMR_contr         2094
## TRIALLowSalinity:SMR_contr     2369
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma    13.27      1.99     8.77    17.08 1.00      442      293
## 
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Conclusions:

• the intercept indicates that the average metabolic change associated with the High Temperature trial (at the mean SMR_contr and MASS) across all fish was 147.67.

• for the same SMR_contr and MASS, the average metabolic change was 107.5 less in the Hypoxia trial and 72.54 less in the Salinity trial than the High Temperature trial.

• in the High Temperature trial, every 1 unit increase in SMR_contr was associated with a decline in metabolic change of round(-1*norin.sum$fixed[4,1], 2) units.

• in the High Temperature trial (and mean SMR_contr), every 1 unit increase in MASS was associated with an increase in metabolic change of round(-1*norin.sum$fixed[5,1], 2) units.

• the rate of change in metabolic change associated with the Hypoxia trial was round(-1*norin.sum$fixed[6,1], 2) units higher than that of the High Temperature trial and round(-1*norin.sum$fixed[7,1], 2) units higher in the Low Salinity trial.

• the variance in intercepts across all fish is 11.85

• the variance in responses between High Temperature and Hypoxia trials across all fish is 6.59.

• the variance in responses between High Temperature and Low Salinity trials across all fish is 28.13.

• the scale of variance between fish within a treatment
  (sigma, 13.27) is substantially higher than that both within the High Temperature trial and between the High Temperature and Hypoxia trials, yet lower than that between High Temperature and Low Salinity trials.

summarise_draws

norin.brm4 %>%
  summarise_draws(
    median,
    ~ HDInterval::hdi(.x),
    rhat,
    ess_bulk,
    ess_tail
  )

## or if you want to exclude some parameters
norin.brm4 %>%
  summarise_draws(
    median,
    ~ HDInterval::hdi(.x),
    rhat,
    ess_bulk,
    ess_tail
  ) %>%
  filter(str_detect(variable, 'prior|^r_|^lp__', negate = TRUE)) 

as_draws_df (posteriors)

norin.brm4 %>% as_draws_df()

norin.brm4 %>%
  as_draws_df() %>%
  summarise_draws(
    median,
    ~ HDInterval::hdi(.x),
    rhat,
    ess_bulk,
    ess_tail
  )

## or if you want to exclude some parameters
norin.brm4 %>%
  as_draws_df() %>%
  summarise_draws(
    median,
    ~ HDInterval::hdi(.x),
    rhat,
    ess_bulk,
    ess_tail
  ) %>%
  filter(str_detect(variable, 'prior|^r_|^lp__', negate = TRUE)) 

tidyMCMC

norin.brm4$fit %>%
    tidyMCMC(estimate.method = 'median',
             conf.int = TRUE,  conf.method = 'HPDinterval',
             rhat = TRUE, ess = TRUE)

Conclusions:

see above

gather_draws

norin.brm4 %>% get_variables()

##   [1] "b_Intercept"                               
##   [2] "b_TRIALHypoxia"                            
##   [3] "b_TRIALLowSalinity"                        
##   [4] "b_SMR_contr"                               
##   [5] "b_MASS"                                    
##   [6] "b_TRIALHypoxia:SMR_contr"                  
##   [7] "b_TRIALLowSalinity:SMR_contr"              
##   [8] "sd_FISHID__Intercept"                      
##   [9] "sd_FISHID__TRIALHypoxia"                   
##  [10] "sd_FISHID__TRIALLowSalinity"               
##  [11] "cor_FISHID__Intercept__TRIALHypoxia"       
##  [12] "cor_FISHID__Intercept__TRIALLowSalinity"   
##  [13] "cor_FISHID__TRIALHypoxia__TRIALLowSalinity"
##  [14] "sigma"                                     
##  [15] "r_FISHID[1,Intercept]"                     
##  [16] "r_FISHID[2,Intercept]"                     
##  [17] "r_FISHID[3,Intercept]"                     
##  [18] "r_FISHID[4,Intercept]"                     
##  [19] "r_FISHID[5,Intercept]"                     
##  [20] "r_FISHID[6,Intercept]"                     
##  [21] "r_FISHID[7,Intercept]"                     
##  [22] "r_FISHID[8,Intercept]"                     
##  [23] "r_FISHID[9,Intercept]"                     
##  [24] "r_FISHID[10,Intercept]"                    
##  [25] "r_FISHID[11,Intercept]"                    
##  [26] "r_FISHID[12,Intercept]"                    
##  [27] "r_FISHID[13,Intercept]"                    
##  [28] "r_FISHID[14,Intercept]"                    
##  [29] "r_FISHID[15,Intercept]"                    
##  [30] "r_FISHID[16,Intercept]"                    
##  [31] "r_FISHID[17,Intercept]"                    
##  [32] "r_FISHID[18,Intercept]"                    
##  [33] "r_FISHID[19,Intercept]"                    
##  [34] "r_FISHID[20,Intercept]"                    
##  [35] "r_FISHID[21,Intercept]"                    
##  [36] "r_FISHID[22,Intercept]"                    
##  [37] "r_FISHID[23,Intercept]"                    
##  [38] "r_FISHID[24,Intercept]"                    
##  [39] "r_FISHID[25,Intercept]"                    
##  [40] "r_FISHID[26,Intercept]"                    
##  [41] "r_FISHID[27,Intercept]"                    
##  [42] "r_FISHID[28,Intercept]"                    
##  [43] "r_FISHID[29,Intercept]"                    
##  [44] "r_FISHID[30,Intercept]"                    
##  [45] "r_FISHID[31,Intercept]"                    
##  [46] "r_FISHID[32,Intercept]"                    
##  [47] "r_FISHID[33,Intercept]"                    
##  [48] "r_FISHID[34,Intercept]"                    
##  [49] "r_FISHID[35,Intercept]"                    
##  [50] "r_FISHID[36,Intercept]"                    
##  [51] "r_FISHID[37,Intercept]"                    
##  [52] "r_FISHID[38,Intercept]"                    
##  [53] "r_FISHID[39,Intercept]"                    
##  [54] "r_FISHID[40,Intercept]"                    
##  [55] "r_FISHID[41,Intercept]"                    
##  [56] "r_FISHID[42,Intercept]"                    
##  [57] "r_FISHID[43,Intercept]"                    
##  [58] "r_FISHID[44,Intercept]"                    
##  [59] "r_FISHID[45,Intercept]"                    
##  [60] "r_FISHID[46,Intercept]"                    
##  [61] "r_FISHID[47,Intercept]"                    
##  [62] "r_FISHID[48,Intercept]"                    
##  [63] "r_FISHID[49,Intercept]"                    
##  [64] "r_FISHID[50,Intercept]"                    
##  [65] "r_FISHID[51,Intercept]"                    
##  [66] "r_FISHID[52,Intercept]"                    
##  [67] "r_FISHID[53,Intercept]"                    
##  [68] "r_FISHID[54,Intercept]"                    
##  [69] "r_FISHID[55,Intercept]"                    
##  [70] "r_FISHID[56,Intercept]"                    
##  [71] "r_FISHID[57,Intercept]"                    
##  [72] "r_FISHID[58,Intercept]"                    
##  [73] "r_FISHID[59,Intercept]"                    
##  [74] "r_FISHID[60,Intercept]"                    
##  [75] "r_FISHID[1,TRIALHypoxia]"                  
##  [76] "r_FISHID[2,TRIALHypoxia]"                  
##  [77] "r_FISHID[3,TRIALHypoxia]"                  
##  [78] "r_FISHID[4,TRIALHypoxia]"                  
##  [79] "r_FISHID[5,TRIALHypoxia]"                  
##  [80] "r_FISHID[6,TRIALHypoxia]"                  
##  [81] "r_FISHID[7,TRIALHypoxia]"                  
##  [82] "r_FISHID[8,TRIALHypoxia]"                  
##  [83] "r_FISHID[9,TRIALHypoxia]"                  
##  [84] "r_FISHID[10,TRIALHypoxia]"                 
##  [85] "r_FISHID[11,TRIALHypoxia]"                 
##  [86] "r_FISHID[12,TRIALHypoxia]"                 
##  [87] "r_FISHID[13,TRIALHypoxia]"                 
##  [88] "r_FISHID[14,TRIALHypoxia]"                 
##  [89] "r_FISHID[15,TRIALHypoxia]"                 
##  [90] "r_FISHID[16,TRIALHypoxia]"                 
##  [91] "r_FISHID[17,TRIALHypoxia]"                 
##  [92] "r_FISHID[18,TRIALHypoxia]"                 
##  [93] "r_FISHID[19,TRIALHypoxia]"                 
##  [94] "r_FISHID[20,TRIALHypoxia]"                 
##  [95] "r_FISHID[21,TRIALHypoxia]"                 
##  [96] "r_FISHID[22,TRIALHypoxia]"                 
##  [97] "r_FISHID[23,TRIALHypoxia]"                 
##  [98] "r_FISHID[24,TRIALHypoxia]"                 
##  [99] "r_FISHID[25,TRIALHypoxia]"                 
## [100] "r_FISHID[26,TRIALHypoxia]"                 
## [101] "r_FISHID[27,TRIALHypoxia]"                 
## [102] "r_FISHID[28,TRIALHypoxia]"                 
## [103] "r_FISHID[29,TRIALHypoxia]"                 
## [104] "r_FISHID[30,TRIALHypoxia]"                 
## [105] "r_FISHID[31,TRIALHypoxia]"                 
## [106] "r_FISHID[32,TRIALHypoxia]"                 
## [107] "r_FISHID[33,TRIALHypoxia]"                 
## [108] "r_FISHID[34,TRIALHypoxia]"                 
## [109] "r_FISHID[35,TRIALHypoxia]"                 
## [110] "r_FISHID[36,TRIALHypoxia]"                 
## [111] "r_FISHID[37,TRIALHypoxia]"                 
## [112] "r_FISHID[38,TRIALHypoxia]"                 
## [113] "r_FISHID[39,TRIALHypoxia]"                 
## [114] "r_FISHID[40,TRIALHypoxia]"                 
## [115] "r_FISHID[41,TRIALHypoxia]"                 
## [116] "r_FISHID[42,TRIALHypoxia]"                 
## [117] "r_FISHID[43,TRIALHypoxia]"                 
## [118] "r_FISHID[44,TRIALHypoxia]"                 
## [119] "r_FISHID[45,TRIALHypoxia]"                 
## [120] "r_FISHID[46,TRIALHypoxia]"                 
## [121] "r_FISHID[47,TRIALHypoxia]"                 
## [122] "r_FISHID[48,TRIALHypoxia]"                 
## [123] "r_FISHID[49,TRIALHypoxia]"                 
## [124] "r_FISHID[50,TRIALHypoxia]"                 
## [125] "r_FISHID[51,TRIALHypoxia]"                 
## [126] "r_FISHID[52,TRIALHypoxia]"                 
## [127] "r_FISHID[53,TRIALHypoxia]"                 
## [128] "r_FISHID[54,TRIALHypoxia]"                 
## [129] "r_FISHID[55,TRIALHypoxia]"                 
## [130] "r_FISHID[56,TRIALHypoxia]"                 
## [131] "r_FISHID[57,TRIALHypoxia]"                 
## [132] "r_FISHID[58,TRIALHypoxia]"                 
## [133] "r_FISHID[59,TRIALHypoxia]"                 
## [134] "r_FISHID[60,TRIALHypoxia]"                 
## [135] "r_FISHID[1,TRIALLowSalinity]"              
## [136] "r_FISHID[2,TRIALLowSalinity]"              
## [137] "r_FISHID[3,TRIALLowSalinity]"              
## [138] "r_FISHID[4,TRIALLowSalinity]"              
## [139] "r_FISHID[5,TRIALLowSalinity]"              
## [140] "r_FISHID[6,TRIALLowSalinity]"              
## [141] "r_FISHID[7,TRIALLowSalinity]"              
## [142] "r_FISHID[8,TRIALLowSalinity]"              
## [143] "r_FISHID[9,TRIALLowSalinity]"              
## [144] "r_FISHID[10,TRIALLowSalinity]"             
## [145] "r_FISHID[11,TRIALLowSalinity]"             
## [146] "r_FISHID[12,TRIALLowSalinity]"             
## [147] "r_FISHID[13,TRIALLowSalinity]"             
## [148] "r_FISHID[14,TRIALLowSalinity]"             
## [149] "r_FISHID[15,TRIALLowSalinity]"             
## [150] "r_FISHID[16,TRIALLowSalinity]"             
## [151] "r_FISHID[17,TRIALLowSalinity]"             
## [152] "r_FISHID[18,TRIALLowSalinity]"             
## [153] "r_FISHID[19,TRIALLowSalinity]"             
## [154] "r_FISHID[20,TRIALLowSalinity]"             
## [155] "r_FISHID[21,TRIALLowSalinity]"             
## [156] "r_FISHID[22,TRIALLowSalinity]"             
## [157] "r_FISHID[23,TRIALLowSalinity]"             
## [158] "r_FISHID[24,TRIALLowSalinity]"             
## [159] "r_FISHID[25,TRIALLowSalinity]"             
## [160] "r_FISHID[26,TRIALLowSalinity]"             
## [161] "r_FISHID[27,TRIALLowSalinity]"             
## [162] "r_FISHID[28,TRIALLowSalinity]"             
## [163] "r_FISHID[29,TRIALLowSalinity]"             
## [164] "r_FISHID[30,TRIALLowSalinity]"             
## [165] "r_FISHID[31,TRIALLowSalinity]"             
## [166] "r_FISHID[32,TRIALLowSalinity]"             
## [167] "r_FISHID[33,TRIALLowSalinity]"             
## [168] "r_FISHID[34,TRIALLowSalinity]"             
## [169] "r_FISHID[35,TRIALLowSalinity]"             
## [170] "r_FISHID[36,TRIALLowSalinity]"             
## [171] "r_FISHID[37,TRIALLowSalinity]"             
## [172] "r_FISHID[38,TRIALLowSalinity]"             
## [173] "r_FISHID[39,TRIALLowSalinity]"             
## [174] "r_FISHID[40,TRIALLowSalinity]"             
## [175] "r_FISHID[41,TRIALLowSalinity]"             
## [176] "r_FISHID[42,TRIALLowSalinity]"             
## [177] "r_FISHID[43,TRIALLowSalinity]"             
## [178] "r_FISHID[44,TRIALLowSalinity]"             
## [179] "r_FISHID[45,TRIALLowSalinity]"             
## [180] "r_FISHID[46,TRIALLowSalinity]"             
## [181] "r_FISHID[47,TRIALLowSalinity]"             
## [182] "r_FISHID[48,TRIALLowSalinity]"             
## [183] "r_FISHID[49,TRIALLowSalinity]"             
## [184] "r_FISHID[50,TRIALLowSalinity]"             
## [185] "r_FISHID[51,TRIALLowSalinity]"             
## [186] "r_FISHID[52,TRIALLowSalinity]"             
## [187] "r_FISHID[53,TRIALLowSalinity]"             
## [188] "r_FISHID[54,TRIALLowSalinity]"             
## [189] "r_FISHID[55,TRIALLowSalinity]"             
## [190] "r_FISHID[56,TRIALLowSalinity]"             
## [191] "r_FISHID[57,TRIALLowSalinity]"             
## [192] "r_FISHID[58,TRIALLowSalinity]"             
## [193] "r_FISHID[59,TRIALLowSalinity]"             
## [194] "r_FISHID[60,TRIALLowSalinity]"             
## [195] "prior_Intercept"                           
## [196] "prior_b_TRIALHypoxia"                      
## [197] "prior_b_TRIALLowSalinity"                  
## [198] "prior_b_SMR_contr"                         
## [199] "prior_b_MASS"                              
## [200] "prior_sigma"                               
## [201] "prior_sd_FISHID"                           
## [202] "prior_cor_FISHID"                          
## [203] "lprior"                                    
## [204] "lp__"                                      
## [205] "accept_stat__"                             
## [206] "treedepth__"                               
## [207] "stepsize__"                                
## [208] "divergent__"                               
## [209] "n_leapfrog__"                              
## [210] "energy__"

norin.draw <- norin.brm4 %>%
    gather_draws(`^b.Intercept$|b_.*|sd_.*|sigma`,  regex=TRUE)
norin.draw

We can then summarise this

norin.draw %>% median_hdci

slab

norin.brm4 %>%
    gather_draws(`b_Intercept|b_.*`, regex=TRUE) %>% 
    ggplot() +
    geom_vline(xintercept=0, linetype='dashed') +
    stat_slab(aes(x = .value, y = .variable,
                  fill = stat(ggdist::cut_cdf_qi(cdf,
                           .width = c(0.5, 0.8, 0.95), 
                           labels = scales::percent_format())
                           )), color='black') + 
    scale_fill_brewer('Interval', direction = -1, na.translate = FALSE) 

norin.brm4 %>% 
    gather_draws(`b_Intercept|b_.*`, regex=TRUE) %>% 
    ggplot() + 
    stat_halfeye(aes(x=.value,  y=.variable)) +
    facet_wrap(~.variable, scales='free')

norin.brm4 %>% 
    gather_draws(`b_Intercept|b_.*`, regex=TRUE) %>% 
    ggplot() + 
    stat_halfeye(aes(x=.value,  y=.variable)) +
    theme_classic()

bayesplot

norin.brm4$fit %>% plot(type='intervals') 

## 'pars' not specified. Showing first 10 parameters by default.

## ci_level: 0.8 (80% intervals)

## outer_level: 0.95 (95% intervals)

half-eye (ggdist)

norin.brm4 %>% 
    gather_draws(`^b_.*`, regex=TRUE) %>% 
    filter(.variable != 'b_Intercept') %>%
    ggplot() + 
    stat_halfeye(aes(x=.value,  y=.variable)) +
    facet_wrap(~.variable, scales='free')

norin.brm4 %>% 
    gather_draws(`^b_.*`, regex=TRUE) %>% 
    filter(.variable != 'b_Intercept') %>%
    ggplot() + 
    stat_halfeye(aes(x=.value,  y=.variable)) +
    geom_vline(xintercept = 0, linetype = 'dashed')

density ridges (ggridges)

norin.brm4 %>% 
    gather_draws(`^b_.*`, regex=TRUE) %>% 
    filter(.variable != 'b_Intercept') %>%
    ggplot() +  
    geom_density_ridges(aes(x=.value, y = .variable), alpha=0.4) +
    geom_vline(xintercept = 0, linetype = 'dashed')

## Picking joint bandwidth of 2.13

##Or in colour
norin.brm4 %>% 
    gather_draws(`^b_.*`, regex=TRUE) %>% 
    filter(.variable != 'b_Intercept') %>%
    ggplot() + 
    geom_density_ridges_gradient(aes(x=exp(.value),
                                     y = .variable,
                                     fill = stat(x)),
                                 alpha=0.4, colour = 'white',
                                 quantile_lines = TRUE,
                                 quantiles = c(0.025, 0.975)) +
    geom_vline(xintercept = 1, linetype = 'dashed') +
    scale_x_continuous(trans = scales::log2_trans()) +
    scale_fill_viridis_c(option = "C") 

## Picking joint bandwidth of 3.07

tidy_draws

norin.brm4 %>% tidy_draws()

spread_draws

norin.brm4 %>% spread_draws(`.*Intercept.*|^b_.*`,  regex=TRUE)

posterior_samples

norin.brm4 %>% posterior_samples() %>% as_tibble()

## Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
## recommended alternatives.

\(R^2\)

norin.brm4 %>%
    bayes_R2(re.form = NA, summary=FALSE) %>%
    median_hdci

norin.brm4 %>%
    bayes_R2(re.form = ~(1|FISHID), summary=FALSE) %>%
    median_hdci

norin.brm4 %>%
    bayes_R2(re.form = ~(TRIAL|FISHID), summary=FALSE) %>%
    median_hdci

ROPE

Region of Practical Equivalence

0.1 * sd(norin$CHANGE)

## [1] 3.38454

norin.brm4 %>% rope(range = c(-3.38, 3.38))

rope(norin.brm4, range = c(-3.38, 3.38)) %>% plot()

Predictions / further analyses

The presence of an interaction suggests that the nature of the relationship between SMR change and the control SMR differs between the three trials (High Temperature, Hypoxia and Low Salinity). Equivalently, the differences between the three trials is not constant over all levels of control SMR.

To further tease these interactions apart, we could either:

• compare the slopes separately for each of the three trials
• compare the effects of the three trials at multiple values of control SMR

We will now explore each of the above, although typically we would only do one or the other (depending on which variable was of more interest).

brms

Compare slopes

norin.brm4 %>% emtrends(~TRIAL, var='SMR_contr') %>% pairs()

##  contrast                      estimate lower.HPD upper.HPD
##  HighTemperature - Hypoxia       -11.06    -18.78     -2.93
##  HighTemperature - LowSalinity    -6.06    -17.78      6.59
##  Hypoxia - LowSalinity             5.11     -9.28     18.50
## 
## Point estimate displayed: median 
## HPD interval probability: 0.95

Conclusions:

• the partial slope associated with control SMR during High Temperature was found to be 11.06 units more negative than during Hypoxia (although not significant). Note, we already had this partial slope in the original summary table.
• the partial slope associated with control SMR during High Temperature was found to be 6.06 units more negative than during Low Salinity (although not significant). Note, we already had this partial slope in the original summary table.
• the partial slope associated with control SMR during Hypoxia was found to be 5.11 units less negative than during Low Salinity (although not significant). Note, we already had this partial slope in the original summary table.

Compare Trials (at different values of control SMR)

We will estimated the pairwise differences between each of the three trials at the minimum, mean and maximum control SMR.

norin.grid <- with(norin,  list(SMR_contr=Hmisc::smean.sdl(SMR_contr)))
norin.grid

## $SMR_contr
##     Mean    Lower    Upper 
## 5.178857 3.926671 6.431044

norin.brm4 %>% emmeans(~TRIAL|SMR_contr,  at=norin.grid) %>% pairs()

## SMR_contr = 3.93:
##  contrast                      estimate lower.HPD upper.HPD
##  HighTemperature - Hypoxia        64.62      53.1    75.570
##  HighTemperature - LowSalinity    49.18      31.2    65.913
##  Hypoxia - LowSalinity           -15.50     -35.0     4.193
## 
## SMR_contr = 5.18:
##  contrast                      estimate lower.HPD upper.HPD
##  HighTemperature - Hypoxia        50.89      44.6    56.069
##  HighTemperature - LowSalinity    41.77      33.3    50.367
##  Hypoxia - LowSalinity            -9.06     -18.6     0.822
## 
## SMR_contr = 6.43:
##  contrast                      estimate lower.HPD upper.HPD
##  HighTemperature - Hypoxia        37.01      25.9    49.858
##  HighTemperature - LowSalinity    34.23      16.8    52.015
##  Hypoxia - LowSalinity            -2.78     -24.4    16.391
## 
## Point estimate displayed: median 
## HPD interval probability: 0.95

norin.brm4 %>%
    emmeans(~TRIAL|SMR_contr,  at=norin.grid) %>%
    pairs() %>%
    gather_emmeans_draws() %>%
    median_hdci()

norin.brm4 %>%
    emmeans(~TRIAL|SMR_contr,  at=norin.grid) %>%
    pairs() %>%
    tidy_draws() %>%
    summarise_draws(median)

Conclusions:

• as control SMR increases, the magnitude of differences between the trials becomes reduced.
• the High Temperature trial resulted in the greatest positive change in SMR.
• the change in SMR was not found to be different between the Hypoxia and Low Salinity trials.

Of course, as with other analyses, it is possible to extract the entire posterior and use it to calculate exceedence probabilities etc.

norin.em <- norin.brm4 %>%
    emmeans(~TRIAL|SMR_contr, at=norin.grid) %>%
    pairs() %>% 
    gather_emmeans_draws() %>%
    mutate(Fit=.value)
norin.em

norin.em %>% group_by(contrast) %>% median_hdci(Fit)

norin.em %>% group_by(contrast, SMR_contr) %>% median_hdci(Fit)

## norin.em %>%
##     group_by(contrast) %>%
##     summarize(P=sum(Fit>0)/n())
norin.em %>%
    group_by(contrast, SMR_contr) %>%
    summarize(P=sum(Fit>0)/n(),
              P2 = 1 - P)

## `summarise()` has grouped output by 'contrast'. You can override using the
## `.groups` argument.

Summary figures

brms

norin.grid <- with(norin, list(SMR_contr = modelr::seq_range(SMR_contr, n = 100)))
newdata <- norin.brm4 %>% emmeans(~SMR_contr|TRIAL, at = norin.grid) %>%
    as.data.frame
head(newdata)

ggplot(data = newdata, aes(y = emmean, x = SMR_contr)) +
  geom_ribbon(aes(ymin = lower.HPD, ymax = upper.HPD, fill = TRIAL), alpha = 0.3) +
  geom_line(aes(, color = TRIAL)) +
  theme_classic() +
  theme(legend.position = c(0.99, 0.99),
        legend.justification = c(1, 1))

## The .fixed values are the predicted values without random effects
obs <- norin.brm4 %>%
  augment() %>%
  mutate(PartialObs=.fitted + .resid)

ggplot(data=newdata, aes(y=emmean, x=SMR_contr)) +
  geom_ribbon(aes(ymin=lower.HPD, ymax=upper.HPD, fill=TRIAL), alpha=0.3) +
  geom_line(aes(, color=TRIAL)) +
  geom_point(data=obs,  aes(y=PartialObs,  color=TRIAL)) +
  ## geom_point(data=norin,  aes(y=CHANGE), color='gray') +
  theme_classic()

g1 <- ggplot(data=newdata, aes(y=emmean, x=SMR_contr)) +
  geom_ribbon(aes(ymin=lower.HPD, ymax=upper.HPD, fill=TRIAL), alpha=0.3) +
  geom_line(aes(, color=TRIAL)) +
  geom_point(data=obs,  aes(y=PartialObs,  color=TRIAL)) +
  theme_classic()
g2 <- ggplot(data=newdata, aes(y=emmean, x=SMR_contr)) +
  geom_ribbon(aes(ymin=lower.HPD, ymax=upper.HPD, fill=TRIAL), alpha=0.3) +
  geom_line(aes(, color=TRIAL)) +
  geom_point(data=norin,  aes(y=CHANGE,  color=TRIAL)) +
  theme_classic()
g1 + g2

References

norin = norin %>% mutate(FISHID=factor(FISHID),
                         TRIAL=factor(TRIAL))

ggplot(norin, aes(y=CHANGE, x=TRIAL)) + geom_boxplot()
ggplot(norin, aes(y=CHANGE, x=SMR_contr, shape=TRIAL, color=TRIAL)) +
    geom_smooth(method='lm') + geom_point()
#ggplot(norin, aes(y=CHANGE, x=MASS, shape=TRIAL, color=TRIAL)) +
#geom_smooth(method='lm') + geom_point()
ggplot(norin, aes(y=CHANGE, x=as.numeric(FISHID), color=TRIAL)) +
    geom_point() + geom_line()

#ggplot(norin, aes(y=MASS, x=TRIAL)) + geom_boxplot()
ggplot(norin, aes(y=CHANGE, x=MASS, color=TRIAL)) + geom_point() + geom_smooth(method='lm')

norin.rstanarm = stan_glmer(CHANGE ~ (1|FISHID)+TRIAL*SMR_contr+MASS, data=norin,
                            prior_PD=TRUE, 
                         iter=5000, warmup=2000, chains=3, thin=5, refresh=0)
prior_summary(norin.rstanarm)

posterior_vs_prior(norin.rstanarm, color_by='vs', group_by=TRUE,
                   facet_args=list(scales='free_y'), pars=c('(Intercept)'))
ggpredict(norin.rstanarm, ~TRIAL*SMR_contr) %>% plot(add.data=TRUE)

norin.rstanarm %>% get_variables()
plot(norin.rstanarm,  'mcmc_trace', regex_pars='^.Intercept|TRIAL|SMR|MASS|[sS]igma')
plot(norin.rstanarm,  'mcmc_acf_bar', regex_pars='^.Intercept|TRIAL|SMR|MASS|[sS]igma')
plot(norin.rstanarm,  'mcmc_rhat_hist', regex_pars='^.Intercept|TRIAL|SMR|MASS|[sS]igma')
plot(norin.rstanarm,  'mcmc_neff_hist', regex_pars='^.Intercept|TRIAL|SMR|MASS|[sS]igma')

#norin.rstan1 = stan_glmer(CHANGE ~ (TRIAL|FISHID)+TRIAL*SMR_contr+MASS, data=norin,
#                          iter=5000, warmup=2000, chains=3, thin=5, refresh=0, cores=3)
norin.rstanarm1 = stan_glmer(CHANGE ~ (SMR_contr|FISHID) + TRIAL*SMR_contr+MASS, data=norin,
                             prior_PD=FALSE, 
                          iter=5000, warmup=2000, chains=3, thin=5, refresh=0, cores=3)
norin.rstanarm1 = update(norin.rstanarm1,  prior_PD=FALSE)



norin.rstanarm2 = stan_glmer(CHANGE ~ (TRIAL*SMR_contr|FISHID) + TRIAL*SMR_contr+MASS, data=norin,
                             prior_PD=FALSE, 
                          iter=5000, warmup=2000, chains=3, thin=5, refresh=0, cores=3)

posterior_vs_prior(norin.rstanarm1, color_by='vs', group_by=TRUE,
                   facet_args=list(scales='free_y'), pars=c('(Intercept)'))

ggpredict(norin.rstanarm1, ~TRIAL*SMR_contr) %>% plot(add.data=TRUE)

norin.rstanarm1 %>% get_variables()
plot(norin.rstanarm1,  'mcmc_trace', regex_pars='^.Intercept|TRIAL|^SMR|MASS|[sS]igma')
plot(norin.rstanarm1,  'mcmc_acf_bar', regex_pars='^.Intercept|TRIAL|^SMR|MASS|[sS]igma')
plot(norin.rstanarm1,  'mcmc_rhat_hist', regex_pars='^.Intercept|TRIAL|^SMR|MASS|[sS]igma')
plot(norin.rstanarm1,  'mcmc_neff_hist', regex_pars='^.Intercept|TRIAL|^SMR|MASS|[sS]igma')


(l.1 <- loo(norin.rstanarm))
(l.2 <- loo(norin.rstanarm1))
loo_compare(l.1,  l.2)


preds <- posterior_predict(norin.rstanarm,  nsamples=250,  summary=FALSE)
norin.resids <- createDHARMa(simulatedResponse = t(preds),
                            observedResponse = norin$CHANGE,
                            fittedPredictedResponse = apply(preds, 2, median))
plot(norin.resids)


g=ggpredict(norin.rstanarm) %>% plot
do.call('grid.arrange', g)

#ggemmeans(norin.rstan, ~TRIAL)

summary(norin.rstanarm)
nms <- norin.rstanarm1 %>% get_variables()
wch <- grep('^.Intercept|TRIAL|^SMR|[sS]igma', nms)
tidyMCMC(norin.rstanarm$stanfit,conf.int=TRUE, conf.method='HPDinterval',
         rhat=TRUE, ess=TRUE, pars=nms[wch], estimate.method='median')

tidyMCMC(norin.rstanarm1$stanfit,conf.int=TRUE, conf.method='HPDinterval',
         rhat=TRUE, ess=TRUE, pars=nms[wch], estimate.method='median')


norin.grid = with(norin, list(SMR_contr=seq(min(SMR_contr),max(SMR_contr), len=100)))
newdata = emmeans(norin.rstanarm, ~TRIAL|SMR_contr, at=norin.grid) %>% as.data.frame
head(newdata)
ggplot(newdata, aes(y=emmean, x=SMR_contr, color=TRIAL)) +
    geom_ribbon(aes(ymin=lower.HPD, ymax=upper.HPD, fill=TRIAL), alpha=0.3,color=NA) +
    geom_line()

norin.grid = with(norin, list(SMR_contr=c(min(SMR_contr),mean(SMR_contr),max(SMR_contr))))

emmeans(norin.rstan, pairwise~TRIAL|SMR_contr, at=norin.grid)

norin.em = emmeans(norin.rstan, pairwise~TRIAL|SMR_contr, at=norin.grid)$contrast %>%
              gather_emmeans_draws() %>%
              mutate(Fit=.value)
norin.em
norin.em %>% group_by(contrast) %>% median_hdci(Fit)
norin.em %>% group_by(contrast, SMR_contr) %>% median_hdci(Fit)
## norin.em %>%
##     group_by(contrast) %>%
##     summarize(P=sum(Fit>0)/n())
norin.em %>%
    group_by(contrast, SMR_contr) %>%
    summarize(P=sum(Fit>0)/n())


bayes_R2(norin.rstanarm, re.form=NA) %>% median_hdi
bayes_R2(norin.rstanarm, re.form=~(1|FISHID)) %>% median_hdi
#bayes_R2(norin.rstan1, re.form=~(SMR_contr|FISHID)) %>% median_hdi

brms

norin = norin %>% mutate(FISHID=factor(FISHID),
                         TRIAL=factor(TRIAL))

ggplot(norin, aes(y=CHANGE, x=TRIAL)) + geom_boxplot()
ggplot(norin, aes(y=CHANGE, x=SMR_contr, shape=TRIAL, color=TRIAL)) +
    geom_smooth(method='lm') + geom_point()
ggplot(norin, aes(y=CHANGE, x=MASS, shape=TRIAL, color=TRIAL)) +
    geom_smooth(method='lm') + geom_point()
ggplot(norin, aes(y=CHANGE, x=as.numeric(FISHID), color=TRIAL)) +
    geom_point() + geom_line()

##ggplot(norin, aes(y=MASS, x=TRIAL)) + geom_boxplot()
##ggplot(norin, aes(y=CHANGE, x=MASS, color=TRIAL)) + geom_point() + geom_smooth(method='lm')

norin %>% group_by(TRIAL) %>%
    summarise(median(CHANGE),
              mad(CHANGE))
priors <- prior(normal(50, 20), class='Intercept') +
    prior(normal(0, 10), class='b') +
    prior(gamma(2,1), class='sigma') +
    prior(gamma(2,1), class='sd')

norin.form <- bf(CHANGE ~ (1|FISHID)+TRIAL*SMR_contr+MASS,
                 family=gaussian)

norin.brm1 = brm(norin.form,
                 data=norin,
                 prior = priors,
                 sample_prior = 'yes',
                 iter=5000, warmup=2000,
                 chains=3, thin=5, refresh=0)

norin.brm1 %>% get_variables()
pars <- norin.brm1 %>% get_variables()
wch <- grepl('^b.Intercept|TRIAL|SMR|MASS|[sS]igma|^sd', pars, perl=TRUE)

stan_trace(norin.brm1$fit, pars=pars[wch])
stan_ac(norin.brm1$fit, pars=pars[wch])
stan_rhat(norin.brm1$fit, pars=pars[wch])
stan_ess(norin.brm1$fit, pars=pars[wch])

##mcmc_plot(norin.brms,  type='trace',
##          regex_pars='^b.*|sigma|^sd')
##mcmc_plot(norin.brms,  type='trace',
##          regex_pars='^b.Intercept|TRIAL|SMR|MASS|[sS]igma|^sd')
##mcmc_plot(norin.brms,  type='acf_bar',
##          regex_pars='^b.Intercept|TRIAL|SMR|MASS|[sS]igma|^sd')
##mcmc_plot(norin.brms,  type='rhat_hist',
##          regex_pars='^b.Intercept|TRIAL|SMR|MASS|[sS]igma|^sd')
##mcmc_plot(norin.brms,  type='neff_hist',
##          regex_pars='^b.Intercept|TRIAL|SMR|MASS|[sS]igma|^sd')

preds <- posterior_predict(norin.brm1,  nsamples=250,  summary=FALSE)
norin.resids <- createDHARMa(simulatedResponse = t(preds),
                            observedResponse = norin$CHANGE,
                            fittedPredictedResponse = apply(preds, 2, median),
                            integerResponse =FALSE)
plot(norin.resids)
#norin.rstan1 = stan_glmer(CHANGE ~ (TRIAL|FISHID)+TRIAL*SMR_contr+MASS, data=norin,
#                          iter=5000, warmup=2000, chains=3, thin=5, refresh=0, cores=3)
norin.form <- bf(CHANGE ~ (TRIAL|FISHID) + TRIAL*SMR_contr+MASS,
                 family=gaussian)
norin.brm2 = brm(norin.form, data=norin,
                 prior = priors,
                 sample_prior = 'yes',
                 iter=5000, warmup=2000,
                 chains=3, thin=5, refresh=0, cores=3,
                 control=list(adapt_delta=0.99))

norin.brm2 %>% get_variables()

pars <- norin.brm2 %>% get_variables()
## wch <- grepl('^b.Intercept|TRIAL|SMR|MASS|[sS]igma|^sd', pars, perl=TRUE)
wch <- grepl('^b_.*|[sS]igma|^sd_.*', pars, perl=TRUE)

stan_trace(norin.brm2$fit, pars=pars[wch])
stan_ac(norin.brm2$fit, pars=pars[wch])
stan_rhat(norin.brm2$fit)#, pars=pars[wch])
stan_ess(norin.brm2$fit)#, pars=pars[wch])
##mcmc_plot(norin.brm2,  type='trace',
##          regex_pars='^b.Intercept|TRIAL|SMR|MASS|[sS]igma|^sd')
##mcmc_plot(norin.brm2,  type='trace',
##          regex_pars='^b.Intercept|TRIAL|SMR|MASS|[sS]igma|^sd')
##mcmc_plot(norin.brm2,  type='acf_bar',
##          regex_pars='^b.Intercept|TRIAL|SMR|MASS|[sS]igma|^sd')
##mcmc_plot(norin.brm2,  type='rhat_hist',
##          regex_pars='^b.Intercept|TRIAL|SMR|MASS|[sS]igma|^sd')
##mcmc_plot(norin.brm2,  type='neff_hist',
##          regex_pars='^b.Intercept|TRIAL|SMR|MASS|[sS]igma|^sd')

(l.1 <- loo(norin.brm1))
(l.2 <- loo(norin.brm2))
loo_compare(l.1,  l.2)


preds <- posterior_predict(norin.brm2,  nsamples=250,  summary=FALSE)
norin.resids <- createDHARMa(simulatedResponse = t(preds),
                            observedResponse = norin$CHANGE,
                            fittedPredictedResponse = apply(preds, 2, median))
plot(norin.resids)

g <- norin.brm2 %>%
    conditional_effects() %>%
    plot(points=TRUE, ask=FALSE)
library(patchwork)
g[[1]] + g[[2]] + g[[3]] + g[[4]]


##g=ggpredict(norin.brms1) %>% plot
##library(patchwork)
##g[[1]] + g[[2]] + g[[3]]

##do.call('grid.arrange', g)

ggemmeans(norin.brm2, ~TRIAL) %>% plot

summary(norin.brm2)

tidyMCMC(norin.brm2$fit,conf.int=TRUE, conf.method='HPDinterval',
         rhat=TRUE, ess=TRUE, estimate.method='median') %>%
  slice(1:11)

pars <- norin.brm2 %>% get_variables()
wch <- grep('^b.Intercept|TRIAL|^b.*SMR|[sS]igma|^sd', pars)
tidyMCMC(norin.brms1$fit,conf.int=TRUE, conf.method='HPDinterval',
         rhat=TRUE, ess=TRUE, pars=pars[wch], estimate.method='median')

bayes_R2(norin.brm2, re.form=NA,  summary=FALSE) %>%
    median_hdci
bayes_R2(norin.brm2, re.form=~(1|FISHID), summary=FALSE) %>%
    median_hdci
bayes_R2(norin.brm2, re.form=~(TRIAL|FISHID), summary=FALSE) %>%
    median_hdci

emmeans(norin.brm2, pairwise~TRIAL)


norin.em <- norin.brm2 %>%
    emmeans(~TRIAL) %>%
    pairs() %>%
    gather_emmeans_draws() %>%
    mutate(Fit=.value)

norin.em %>%
  group_by(contrast) %>%
  median_hdi()

norin.em %>%
    ggplot() +
    geom_vline(xintercept=0, linetype='dashed') +
    stat_slab(aes(x=.value, y=contrast,
                  fill = stat(ggdist::cut_cdf_qi(cdf,
                            .width = c(0.5, 0.8, 0.95), 
                            labels = scales::percent_format())
                            )), color='black') +
    scale_fill_brewer('Interval', direction = -1, na.translate = FALSE) +
    theme_bw()

norin.em %>%
    group_by(contrast) %>%
  summarize(P=sum(Fit>0)/n())


norin.grid <- with(norin, list(SMR_contr=c(min(SMR_contr),
                                           mean(SMR_contr),
                                           max(SMR_contr))))

norin.em <- norin.brm2 %>%
    emmeans(~TRIAL|SMR_contr, at=norin.grid) %>%
    pairs() %>%
    gather_emmeans_draws()

norin.em %>% head
norin.em %>%
    group_by(contrast, SMR_contr) %>%
    median_hdi()

norin.em %>%
    group_by(contrast, SMR_contr) %>%
    summarize(P=sum(.value>0)/n())

norin.grid <- with(norin, list(SMR_contr=modelr::seq_range(SMR_contr, n=100)))
newdata <- norin.brm2 %>%
    emmeans(~SMR_contr|TRIAL, at=norin.grid) %>%
    as.data.frame
head(newdata)
partial.obs <- norin %>%
    mutate(Pred = predict(norin.brm2, re.form = NA, summary=TRUE)[,'Estimate'],
           Resid = resid(norin.brm2)[,'Estimate'],
           Obs = Pred + Resid)
ggplot(newdata, aes(y=emmean, x=SMR_contr, color=TRIAL)) +
    geom_point(data=partial.obs, aes(y=Obs)) +
    ##geom_point(data=partial.obs, aes(y=CHANGE), shape=2) +
    geom_ribbon(aes(ymin=lower.HPD, ymax=upper.HPD, fill=TRIAL), alpha=0.3,color=NA) +
    geom_line()

Norin, Tommy, Hans Malte, and Timothy D. Clark. 2015. “Differential Plasticity of Metabolic Rate Phenotypes in a Tropical Fish Facing Environmental Change.” Functional Ecology 30 (3): 369–78. https://doi.org/10.1111/1365-2435.12503.